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.

agrif_oce_update.F90 in NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/NST/agrif_oce_update.F90 @ 13680

Last change on this file since 13680 was 13680, checked in by jchanut, 4 years ago
#2385, effectively add VORTEX ssh initialization and key_qco. Save first guess (prior update) r3x arrays: Restored correct value of surface w.
Property svn:keywords set to `Id`
File size: 58.8 KB

Line
1	#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */
2	#undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */
3	#undef VOL_REFLUX /* VOLUME REFLUXING*/
4
5	MODULE agrif_oce_update
6	!!======================================================================
7	!! * MODULE agrif_oce_interp *
8	!! AGRIF: update package for the ocean dynamics (OPA)
9	!!======================================================================
10	!! History : 2.0 ! 2002-06 (L. Debreu) Original code
11	!! 3.2 ! 2009-04 (R. Benshila)
12	!! 3.6 ! 2014-09 (R. Benshila)
13	!!----------------------------------------------------------------------
14	#if defined key_agrif
15	!!----------------------------------------------------------------------
16	!! 'key_agrif' AGRIF zoom
17	!!----------------------------------------------------------------------
18	USE par_oce
19	USE oce
20	USE dom_oce
21	USE zdf_oce ! vertical physics: ocean variables
22	USE agrif_oce
23	!
24	USE in_out_manager ! I/O manager
25	USE lib_mpp ! MPP library
26	USE domvvl ! Need interpolation routines
27	USE vremap ! Vertical remapping
28	USE lbclnk
29	#if defined key_qco
30	USE domqco
31	#endif
32	IMPLICIT NONE
33	PRIVATE
34
35	PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh
36	PUBLIC Update_Scales
37
38	!! * Substitutions
39	# include "domzgr_substitute.h90"
40	!!----------------------------------------------------------------------
41	!! NEMO/NST 4.0 , NEMO Consortium (2018)
42	!! $Id$
43	!! Software governed by the CeCILL license (see ./LICENSE)
44	!!----------------------------------------------------------------------
45	CONTAINS
46
47	SUBROUTINE Agrif_Update_Tra( )
48	!!----------------------------------------------------------------------
49	!! * ROUTINE Agrif_Update_Tra *
50	!!----------------------------------------------------------------------
51	!
52	IF (Agrif_Root()) RETURN
53	!
54	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update tracers from grid Number',Agrif_Fixed()
55
56	#if defined key_vertical
57	! Effect of this has to be carrefully checked
58	! depending on what the nesting tools ensure for
59	! volume conservation:
60	Agrif_UseSpecialValueInUpdate = .FALSE.
61	#else
62	Agrif_UseSpecialValueInUpdate = .TRUE.
63	#endif
64	Agrif_SpecialValueFineGrid = 0._wp
65	!
66	# if ! defined DECAL_FEEDBACK
67	CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
68	! near boundary update:
69	! CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
70	# else
71	CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
72	! near boundary update:
73	! CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
74	# endif
75	!
76	Agrif_UseSpecialValueInUpdate = .FALSE.
77	!
78	!
79	END SUBROUTINE Agrif_Update_Tra
80
81	SUBROUTINE Agrif_Update_Dyn( )
82	!!----------------------------------------------------------------------
83	!! * ROUTINE Agrif_Update_Dyn *
84	!!----------------------------------------------------------------------
85	!
86	IF (Agrif_Root()) RETURN
87	!
88	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update momentum from grid Number',Agrif_Fixed()
89
90	Agrif_UseSpecialValueInUpdate = .FALSE.
91	Agrif_SpecialValueFineGrid = 0._wp
92
93	use_sign_north = .TRUE.
94	sign_north = -1._wp
95
96	!
97	# if ! defined DECAL_FEEDBACK
98	CALL Agrif_Update_Variable(un_update_id,procname = updateU)
99	CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
100	! near boundary update:
101	! CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
102	! CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
103	# else
104	CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
105	CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
106	! near boundary update:
107	! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
108	! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
109	# endif
110
111	# if ! defined DECAL_FEEDBACK_2D
112	CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
113	CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)
114	# else
115	CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
116	CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)
117	# endif
118	!
119	# if ! defined DECAL_FEEDBACK_2D
120	! Account for updated thicknesses at boundary edges
121	IF (.NOT.ln_linssh) THEN
122	! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
123	! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
124	ENDIF
125	# endif
126	!
127	IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
128	! Update time integrated transports
129	# if ! defined DECAL_FEEDBACK_2D
130	CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
131	CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
132	# else
133	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
134	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
135	# endif
136	END IF
137	!
138	use_sign_north = .FALSE.
139	!
140	END SUBROUTINE Agrif_Update_Dyn
141
142	SUBROUTINE Agrif_Update_ssh( )
143	!!---------------------------------------------
144	!! * ROUTINE Agrif_Update_ssh *
145	!!---------------------------------------------
146	!
147	IF (Agrif_Root()) RETURN
148	!
149	Agrif_UseSpecialValueInUpdate = .TRUE.
150	Agrif_SpecialValueFineGrid = 0._wp
151	# if ! defined DECAL_FEEDBACK_2D
152	CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
153	# else
154	CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
155	# endif
156	!
157	Agrif_UseSpecialValueInUpdate = .FALSE.
158	!
159	# if defined VOL_REFLUX
160	IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
161	use_sign_north = .TRUE.
162	sign_north = -1._wp
163	! Refluxing on ssh:
164	# if defined DECAL_FEEDBACK_2D
165	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
166	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
167	# else
168	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu)
169	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv)
170	# endif
171	use_sign_north = .FALSE.
172	END IF
173	# endif
174	!
175	END SUBROUTINE Agrif_Update_ssh
176
177
178	SUBROUTINE Agrif_Update_Tke( )
179	!!---------------------------------------------
180	!! * ROUTINE Agrif_Update_Tke *
181	!!---------------------------------------------
182	!!
183	!
184	IF (Agrif_Root()) RETURN
185	!
186	Agrif_UseSpecialValueInUpdate = .TRUE.
187	Agrif_SpecialValueFineGrid = 0.
188
189	CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN )
190	CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
191	CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
192
193	Agrif_UseSpecialValueInUpdate = .FALSE.
194
195	END SUBROUTINE Agrif_Update_Tke
196
197	SUBROUTINE Agrif_Update_vvl( )
198	!!---------------------------------------------
199	!! * ROUTINE Agrif_Update_vvl *
200	!!---------------------------------------------
201	!
202	IF (Agrif_Root()) RETURN
203	!
204	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
205	!
206	#if ! defined key_qco
207	Agrif_UseSpecialValueInUpdate = .TRUE.
208	Agrif_SpecialValueFineGrid = 0.
209	!
210	! No interface separation here, update vertical grid at T points
211	! everywhere over the overlapping regions (one account for refluxing in that case):
212	CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)
213	!
214	Agrif_UseSpecialValueInUpdate = .FALSE.
215	!
216	CALL Agrif_ChildGrid_To_ParentGrid()
217	CALL dom_vvl_update_UVF
218	CALL Agrif_ParentGrid_To_ChildGrid()
219	#else
220	CALL Agrif_ChildGrid_To_ParentGrid()
221	r3t(:,:,Krhs_a) = r3t(:,:,Kmm_a)
222	r3u(:,:,Krhs_a) = r3u(:,:,Kmm_a)
223	r3v(:,:,Krhs_a) = r3v(:,:,Kmm_a)
224	CALL dom_qco_zgr( Kbb_a, Kmm_a, Krhs_a )
225	CALL Agrif_ParentGrid_To_ChildGrid()
226	#endif
227	!
228	END SUBROUTINE Agrif_Update_vvl
229
230	#if ! defined key_qco
231	SUBROUTINE dom_vvl_update_UVF
232	!!---------------------------------------------
233	!! * ROUTINE dom_vvl_update_UVF *
234	!!---------------------------------------------
235	!!
236	INTEGER :: jk
237	REAL(wp):: zcoef
238	!!---------------------------------------------
239	IF (lwp.AND.lk_agrif_debug) Write(,) 'Finalize e3 on grid Number', &
240	& Agrif_Fixed(), 'Step', Agrif_Nb_Step()
241
242	! Save "old" scale factor (prior update) for subsequent asselin correction
243	! of prognostic variables
244	! -----------------------
245	!
246	e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a)
247	e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a)
248	! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a)
249	! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a)
250	hu(:,:,Krhs_a) = hu(:,:,Kmm_a)
251	hv(:,:,Krhs_a) = hv(:,:,Kmm_a)
252
253	! 1) NOW fields
254	!--------------
255
256	! Vertical scale factor interpolations
257	! ------------------------------------
258	CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) , 'U' )
259	CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) , 'V' )
260	CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) , 'F' )
261
262	CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' )
263	CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' )
264
265	! Update total depths:
266	! --------------------
267	hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points
268	hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points
269	DO jk = 1, jpkm1
270	hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)
271	hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk)
272	END DO
273	! ! Inverse of the local depth
274	r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) )
275	r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) )
276
277
278	! 2) BEFORE fields:
279	!------------------
280	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
281	!
282	! Vertical scale factor interpolations
283	! ------------------------------------
284	CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a), 'U' )
285	CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a), 'V' )
286
287	CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' )
288	CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' )
289
290	! Update total depths:
291	! --------------------
292	hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points
293	hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points
294	DO jk = 1, jpkm1
295	hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk)
296	hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk)
297	END DO
298	! ! Inverse of the local depth
299	r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) )
300	r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) )
301	ENDIF
302	!
303	END SUBROUTINE dom_vvl_update_UVF
304	#endif
305
306	#if defined key_vertical
307
308	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
309	!!----------------------------------------------------------------------
310	!! * ROUTINE updateT *
311	!!---------------------------------------------
312	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
313	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
314	LOGICAL, INTENT(in) :: before
315	!!
316	INTEGER :: ji,jj,jk,jn
317	INTEGER :: N_in, N_out
318	REAL(wp) :: ztb, ztnu, ztno
319	REAL(wp) :: h_in(k1:k2)
320	REAL(wp) :: h_out(1:jpk)
321	REAL(wp) :: tabin(k1:k2,1:jpts)
322	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child
323	!!---------------------------------------------
324	!
325	IF (before) THEN
326	!jc_alt
327	! AGRIF_SpecialValue = -999._wp
328	DO jn = n1,n2-1
329	DO jk=k1,k2
330	DO jj=j1,j2
331	DO ji=i1,i2
332	!jc_alt
333	! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
334	! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp
335	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
336	END DO
337	END DO
338	END DO
339	END DO
340	DO jk=k1,k2
341	DO jj=j1,j2
342	DO ji=i1,i2
343	!jc_alt
344	! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
345	! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp
346	tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
347	END DO
348	END DO
349	END DO
350	ELSE
351	tabres_child(:,:,:,:) = 0._wp
352	AGRIF_SpecialValue = 0._wp
353	DO jj=j1,j2
354	DO ji=i1,i2
355	N_in = 0
356	DO jk=k1,k2 !k2 = jpk of child grid
357	! jc_alt
358	! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT
359	IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT
360	N_in = N_in + 1
361	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
362	h_in(N_in) = tabres(ji,jj,jk,n2)
363	ENDDO
364	N_out = 0
365	DO jk=1,jpk ! jpk of parent grid
366	IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
367	N_out = N_out + 1
368	h_out(N_out) = e3t(ji,jj,jk,Kmm_a)
369	ENDDO
370	IF (N_in*N_out > 0) THEN !Remove this?
371	CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
372	ENDIF
373	ENDDO
374	ENDDO
375
376	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
377	! Add asselin part
378	DO jn = 1,jpts
379	DO jk = 1, jpkm1
380	DO jj = j1, j2
381	DO ji = i1, i2
382	IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
383	ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
384	ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
385	ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
386	ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &
387	& * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
388	ENDIF
389	END DO
390	END DO
391	END DO
392	END DO
393	ENDIF
394	DO jn = 1,jpts
395	DO jk = 1, jpkm1
396	DO jj = j1, j2
397	DO ji = i1, i2
398	IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
399	ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
400	END IF
401	END DO
402	END DO
403	END DO
404	END DO
405	!
406	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
407	ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a)
408	ENDIF
409	ENDIF
410	!
411	END SUBROUTINE updateTS
412
413	# else
414
415	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
416	!!---------------------------------------------
417	!! * ROUTINE updateT *
418	!!---------------------------------------------
419	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
420	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
421	LOGICAL, INTENT(in) :: before
422	!!
423	INTEGER :: ji,jj,jk,jn
424	REAL(wp) :: ztb, ztnu, ztno
425	!!---------------------------------------------
426	!
427	IF (before) THEN
428	DO jn = 1,jpts
429	DO jk=k1,k2
430	DO jj=j1,j2
431	DO ji=i1,i2
432	!> jc tmp
433	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
434	! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
435	!< jc tmp
436	END DO
437	END DO
438	END DO
439	END DO
440	ELSE
441	!> jc tmp
442	DO jn = 1,jpts
443	tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
444	& * tmask(i1:i2,j1:j2,k1:k2)
445	ENDDO
446	!< jc tmp
447	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
448	! Add asselin part
449	DO jn = 1,jpts
450	DO jk = k1, k2
451	DO jj = j1, j2
452	DO ji = i1, i2
453	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
454	ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
455	ztnu = tabres(ji,jj,jk,jn)
456	ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
457	ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &
458	& * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
459	ENDIF
460	END DO
461	END DO
462	END DO
463	END DO
464	ENDIF
465	DO jn = 1,jpts
466	DO jk=k1,k2
467	DO jj=j1,j2
468	DO ji=i1,i2
469	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
470	ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
471	END IF
472	END DO
473	END DO
474	END DO
475	END DO
476	!
477	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
478	ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)
479	ENDIF
480	!
481	ENDIF
482	!
483	END SUBROUTINE updateTS
484
485	# endif
486
487	# if defined key_vertical
488
489	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
490	!!---------------------------------------------
491	!! * ROUTINE updateu *
492	!!---------------------------------------------
493	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
494	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
495	LOGICAL , INTENT(in ) :: before
496	!
497	INTEGER :: ji, jj, jk
498	REAL(wp):: zrhoy, zub, zunu, zuno
499	! VERTICAL REFINEMENT BEGIN
500	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
501	REAL(wp) :: h_in(k1:k2)
502	REAL(wp) :: h_out(1:jpk)
503	INTEGER :: N_in, N_out
504	REAL(wp) :: h_diff, excess, thick
505	REAL(wp) :: tabin(k1:k2)
506	! VERTICAL REFINEMENT END
507	!!---------------------------------------------
508	!
509	IF( before ) THEN
510	zrhoy = Agrif_Rhoy()
511	!jc_alt
512	! AGRIF_SpecialValue = -999._wp
513	DO jk=k1,k2
514	DO jj=j1,j2
515	DO ji=i1,i2
516	!jc_alt
517	! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) &
518	! & + (umask(ji,jj,jk)-1._wp)*999._wp
519	tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)
520	!jc_alt
521	! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) &
522	! & + (umask(ji,jj,jk)-1._wp)*999._wp
523	tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
524	END DO
525	END DO
526	END DO
527	ELSE
528	tabres_child(:,:,:) = 0.
529	AGRIF_SpecialValue = 0._wp
530	DO jj=j1,j2
531	DO ji=i1,i2
532	N_in = 0
533	h_in(:) = 0._wp
534	tabin(:) = 0._wp
535	DO jk=k1,k2 !k2=jpk of child grid
536	!jc_alt
537	! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT
538	IF( tabres(ji,jj,jk,2) == 0.) EXIT
539	N_in = N_in + 1
540	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
541	h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
542	ENDDO
543	N_out = 0
544	DO jk=1,jpk
545	IF (umask(ji,jj,jk) == 0) EXIT
546	N_out = N_out + 1
547	h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
548	ENDDO
549	IF (N_in * N_out > 0) THEN
550	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
551	excess = 0._wp
552	IF (h_diff < -1.e-4) THEN
553	!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
554	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
555	DO jk=N_in,1,-1
556	thick = MIN(-1*h_diff, h_in(jk))
557	excess = excess + tabin(jk)thicke2u(ji,jj)
558	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
559	h_diff = h_diff + thick
560	IF ( h_diff == 0) THEN
561	N_in = jk
562	h_in(jk) = h_in(jk) - thick
563	EXIT
564	ENDIF
565	ENDDO
566	ENDIF
567	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
568	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
569	ENDIF
570	ENDDO
571	ENDDO
572	!
573	DO jk=1,jpk
574	DO jj=j1,j2
575	DO ji=i1,i2
576	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
577	zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
578	zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
579	zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a)
580	uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &
581	& * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
582	ENDIF
583	!
584	uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
585	END DO
586	END DO
587	END DO
588	!
589	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
590	uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a)
591	ENDIF
592	!
593	ENDIF
594	!
595	END SUBROUTINE updateu
596
597	#else
598
599	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
600	!!---------------------------------------------
601	!! * ROUTINE updateu *
602	!!---------------------------------------------
603	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
604	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
605	LOGICAL , INTENT(in ) :: before
606	!
607	INTEGER :: ji, jj, jk
608	REAL(wp) :: zrhoy, zub, zunu, zuno
609	!!---------------------------------------------
610	!
611	IF( before ) THEN
612	zrhoy = Agrif_Rhoy()
613	DO jk = k1, k2
614	tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a)
615	END DO
616	ELSE
617	DO jk=k1,k2
618	DO jj=j1,j2
619	DO ji=i1,i2
620	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)
621	!
622	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
623	zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
624	zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
625	zunu = tabres(ji,jj,jk,1)
626	uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &
627	& * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
628	ENDIF
629	!
630	uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a)
631	END DO
632	END DO
633	END DO
634	!
635	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
636	uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a)
637	ENDIF
638	!
639	ENDIF
640	!
641	END SUBROUTINE updateu
642
643	# endif
644
645	SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
646	!!---------------------------------------------
647	!! * ROUTINE correct_u_bdy *
648	!!---------------------------------------------
649	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
650	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
651	LOGICAL , INTENT(in ) :: before
652	INTEGER , INTENT(in) :: nb, ndir
653	!!
654	LOGICAL :: western_side, eastern_side
655	!
656	INTEGER :: jj, jk
657	REAL(wp) :: zcor
658	!!---------------------------------------------
659	!
660	IF( .NOT.before ) THEN
661	!
662	western_side = (nb == 1).AND.(ndir == 1)
663	eastern_side = (nb == 1).AND.(ndir == 2)
664	!
665	IF (western_side) THEN
666	DO jj=j1,j2
667	zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a)
668	uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor
669	DO jk=1,jpkm1
670	uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk)
671	END DO
672	END DO
673	ENDIF
674	!
675	IF (eastern_side) THEN
676	DO jj=j1,j2
677	zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a)
678	uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor
679	DO jk=1,jpkm1
680	uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk)
681	END DO
682	END DO
683	ENDIF
684	!
685	ENDIF
686	!
687	END SUBROUTINE correct_u_bdy
688
689	# if defined key_vertical
690
691	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
692	!!---------------------------------------------
693	!! * ROUTINE updatev *
694	!!---------------------------------------------
695	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
696	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
697	LOGICAL , INTENT(in ) :: before
698	!
699	INTEGER :: ji, jj, jk
700	REAL(wp) :: zrhox, zvb, zvnu, zvno
701	! VERTICAL REFINEMENT BEGIN
702	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
703	REAL(wp) :: h_in(k1:k2)
704	REAL(wp) :: h_out(1:jpk)
705	INTEGER :: N_in, N_out
706	REAL(wp) :: h_diff, excess, thick
707	REAL(wp) :: tabin(k1:k2)
708	! VERTICAL REFINEMENT END
709	!!---------------------------------------------
710	!
711	IF( before ) THEN
712	zrhox = Agrif_Rhox()
713	!jc_alt
714	! AGRIF_SpecialValue = -999._wp
715	DO jk=k1,k2
716	DO jj=j1,j2
717	DO ji=i1,i2
718	!jc_alt
719	! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) &
720	! & + (vmask(ji,jj,jk)-1._wp) * 999._wp
721	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a)
722	!jc_alt
723	! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) &
724	! & + (vmask(ji,jj,jk)-1._wp) * 999._wp
725	tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk)
726	END DO
727	END DO
728	END DO
729	ELSE
730	tabres_child(:,:,:) = 0.
731	AGRIF_SpecialValue = 0._wp
732	DO jj=j1,j2
733	DO ji=i1,i2
734	N_in = 0
735	DO jk=k1,k2
736	!jc_alt
737	! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT
738	IF (tabres(ji,jj,jk,2) == 0) EXIT
739	N_in = N_in + 1
740	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
741	h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
742	ENDDO
743	N_out = 0
744	DO jk=1,jpk
745	IF (vmask(ji,jj,jk) == 0) EXIT
746	N_out = N_out + 1
747	h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
748	ENDDO
749	IF (N_in * N_out > 0) THEN
750	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
751	excess = 0._wp
752	IF (h_diff < -1.e-4) then
753	!Even if bathy at T points match it's possible for the V points to be deeper in the child grid.
754	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
755	DO jk=N_in,1,-1
756	thick = MIN(-1*h_diff, h_in(jk))
757	excess = excess + tabin(jk)thicke2u(ji,jj)
758	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
759	h_diff = h_diff + thick
760	IF ( h_diff == 0) THEN
761	N_in = jk
762	h_in(jk) = h_in(jk) - thick
763	EXIT
764	ENDIF
765	ENDDO
766	ENDIF
767	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
768	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
769	ENDIF
770	ENDDO
771	ENDDO
772	!
773	DO jk=1,jpkm1
774	DO jj=j1,j2
775	DO ji=i1,i2
776	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
777	zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
778	zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
779	zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a)
780	vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &
781	& * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
782	ENDIF
783	!
784	vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
785	END DO
786	END DO
787	END DO
788	!
789	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
790	vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a)
791	ENDIF
792	!
793	ENDIF
794	!
795	END SUBROUTINE updatev
796
797	# else
798
799	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
800	!!---------------------------------------------
801	!! * ROUTINE updatev *
802	!!---------------------------------------------
803	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
804	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
805	LOGICAL , INTENT(in ) :: before
806	!
807	INTEGER :: ji, jj, jk
808	REAL(wp) :: zrhox, zvb, zvnu, zvno
809	!!---------------------------------------------
810	!
811	IF (before) THEN
812	zrhox = Agrif_Rhox()
813	DO jk=k1,k2
814	DO jj=j1,j2
815	DO ji=i1,i2
816	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
817	END DO
818	END DO
819	END DO
820	ELSE
821	DO jk=k1,k2
822	DO jj=j1,j2
823	DO ji=i1,i2
824	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
825	!
826	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
827	zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
828	zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
829	zvnu = tabres(ji,jj,jk,1)
830	vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &
831	& * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
832	ENDIF
833	!
834	vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a)
835	END DO
836	END DO
837	END DO
838	!
839	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
840	vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a)
841	ENDIF
842	!
843	ENDIF
844	!
845	END SUBROUTINE updatev
846
847	# endif
848
849	SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
850	!!---------------------------------------------
851	!! * ROUTINE correct_v_bdy *
852	!!---------------------------------------------
853	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
854	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
855	LOGICAL , INTENT(in ) :: before
856	INTEGER , INTENT(in) :: nb, ndir
857	!!
858	LOGICAL :: southern_side, northern_side
859	!
860	INTEGER :: ji, jk
861	REAL(wp) :: zcor
862	!!---------------------------------------------
863	!
864	IF( .NOT.before ) THEN
865	!
866	southern_side = (nb == 2).AND.(ndir == 1)
867	northern_side = (nb == 2).AND.(ndir == 2)
868	!
869	IF (southern_side) THEN
870	DO ji=i1,i2
871	zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a)
872	vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor
873	DO jk=1,jpkm1
874	vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk)
875	END DO
876	END DO
877	ENDIF
878	!
879	IF (northern_side) THEN
880	DO ji=i1,i2
881	zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a)
882	vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor
883	DO jk=1,jpkm1
884	vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk)
885	END DO
886	END DO
887	ENDIF
888	!
889	ENDIF
890	!
891	END SUBROUTINE correct_v_bdy
892
893
894	SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
895	!!----------------------------------------------------------------------
896	!! * ROUTINE updateu2d *
897	!!----------------------------------------------------------------------
898	INTEGER , INTENT(in ) :: i1, i2, j1, j2
899	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
900	LOGICAL , INTENT(in ) :: before
901	!!
902	INTEGER :: ji, jj, jk
903	REAL(wp) :: zrhoy
904	REAL(wp) :: zcorr
905	!!---------------------------------------------
906	!
907	IF( before ) THEN
908	zrhoy = Agrif_Rhoy()
909	DO jj=j1,j2
910	DO ji=i1,i2
911	tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj)
912	END DO
913	END DO
914	ELSE
915	DO jj=j1,j2
916	DO ji=i1,i2
917	tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj)
918	!
919	! Update "now" 3d velocities:
920	spgu(ji,jj) = 0._wp
921	DO jk=1,jpkm1
922	spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
923	END DO
924	!
925	zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)
926	DO jk=1,jpkm1
927	uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)
928	END DO
929	!
930	! Update barotropic velocities:
931	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
932	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
933	zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a)
934	uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1)
935	END IF
936	ENDIF
937	uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1)
938	!
939	! Correct "before" velocities to hold correct bt component:
940	spgu(ji,jj) = 0.e0
941	DO jk=1,jpkm1
942	spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
943	END DO
944	!
945	zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a)
946	DO jk=1,jpkm1
947	uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)
948	END DO
949	!
950	END DO
951	END DO
952	!
953	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
954	uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a)
955	ENDIF
956	ENDIF
957	!
958	END SUBROUTINE updateu2d
959
960
961	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
962	!!----------------------------------------------------------------------
963	!! * ROUTINE updatev2d *
964	!!----------------------------------------------------------------------
965	INTEGER , INTENT(in ) :: i1, i2, j1, j2
966	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
967	LOGICAL , INTENT(in ) :: before
968	!
969	INTEGER :: ji, jj, jk
970	REAL(wp) :: zrhox, zcorr
971	!!----------------------------------------------------------------------
972	!
973	IF( before ) THEN
974	zrhox = Agrif_Rhox()
975	DO jj=j1,j2
976	DO ji=i1,i2
977	tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj)
978	END DO
979	END DO
980	ELSE
981	DO jj=j1,j2
982	DO ji=i1,i2
983	tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj)
984	!
985	! Update "now" 3d velocities:
986	spgv(ji,jj) = 0.e0
987	DO jk=1,jpkm1
988	spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
989	END DO
990	!
991	zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)
992	DO jk=1,jpkm1
993	vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)
994	END DO
995	!
996	! Update barotropic velocities:
997	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
998	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
999	zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a)
1000	vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1)
1001	END IF
1002	ENDIF
1003	vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1)
1004	!
1005	! Correct "before" velocities to hold correct bt component:
1006	spgv(ji,jj) = 0.e0
1007	DO jk=1,jpkm1
1008	spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
1009	END DO
1010	!
1011	zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a)
1012	DO jk=1,jpkm1
1013	vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)
1014	END DO
1015	!
1016	END DO
1017	END DO
1018	!
1019	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1020	vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a)
1021	ENDIF
1022	!
1023	ENDIF
1024	!
1025	END SUBROUTINE updatev2d
1026
1027
1028	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
1029	!!----------------------------------------------------------------------
1030	!! * ROUTINE updateSSH *
1031	!!----------------------------------------------------------------------
1032	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1033	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1034	LOGICAL , INTENT(in ) :: before
1035	!!
1036	INTEGER :: ji, jj
1037	!!----------------------------------------------------------------------
1038	!
1039	IF( before ) THEN
1040	DO jj=j1,j2
1041	DO ji=i1,i2
1042	tabres(ji,jj) = ssh(ji,jj,Kmm_a)
1043	END DO
1044	END DO
1045	ELSE
1046	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
1047	DO jj=j1,j2
1048	DO ji=i1,i2
1049	ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) &
1050	& + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
1051	END DO
1052	END DO
1053	ENDIF
1054	!
1055	DO jj=j1,j2
1056	DO ji=i1,i2
1057	ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1)
1058	END DO
1059	END DO
1060	!
1061	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1062	ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a)
1063	ENDIF
1064	!
1065
1066	ENDIF
1067	!
1068	END SUBROUTINE updateSSH
1069
1070
1071	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1072	!!----------------------------------------------------------------------
1073	!! * ROUTINE updateub2b *
1074	!!----------------------------------------------------------------------
1075	INTEGER , INTENT(in) :: i1, i2, j1, j2
1076	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1077	LOGICAL , INTENT(in) :: before
1078	!!
1079	INTEGER :: ji, jj
1080	REAL(wp) :: zrhoy, za1, zcor
1081	!!---------------------------------------------
1082	!
1083	IF (before) THEN
1084	zrhoy = Agrif_Rhoy()
1085	DO jj=j1,j2
1086	DO ji=i1,i2
1087	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1088	END DO
1089	END DO
1090	tabres = zrhoy * tabres
1091	ELSE
1092	!
1093	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1094	!
1095	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1096	DO jj=j1,j2
1097	DO ji=i1,i2
1098	zcor=tabres(ji,jj) - ub2_b(ji,jj)
1099	! Update time integrated fluxes also in case of multiply nested grids:
1100	ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor
1101	! Update corrective fluxes:
1102	un_bf(ji,jj) = un_bf(ji,jj) + zcor
1103	! Update half step back fluxes:
1104	ub2_b(ji,jj) = tabres(ji,jj)
1105	END DO
1106	END DO
1107	ENDIF
1108	!
1109	END SUBROUTINE updateub2b
1110
1111	SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1112	!!---------------------------------------------
1113	!! * ROUTINE reflux_sshu *
1114	!!---------------------------------------------
1115	INTEGER, INTENT(in) :: i1, i2, j1, j2
1116	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1117	LOGICAL, INTENT(in) :: before
1118	INTEGER, INTENT(in) :: nb, ndir
1119	!!
1120	LOGICAL :: western_side, eastern_side
1121	INTEGER :: ji, jj
1122	REAL(wp) :: zrhoy, za1, zcor
1123	!!---------------------------------------------
1124	!
1125	IF (before) THEN
1126	zrhoy = Agrif_Rhoy()
1127	DO jj=j1,j2
1128	DO ji=i1,i2
1129	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1130	END DO
1131	END DO
1132	tabres = zrhoy * tabres
1133	ELSE
1134	!
1135	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1136	!
1137	western_side = (nb == 1).AND.(ndir == 1)
1138	eastern_side = (nb == 1).AND.(ndir == 2)
1139	!
1140	IF (western_side) THEN
1141	DO jj=j1,j2
1142	zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
1143	ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor
1144	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + rn_atfp * zcor
1145	END DO
1146	ENDIF
1147	IF (eastern_side) THEN
1148	DO jj=j1,j2
1149	zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1150	ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
1151	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor
1152	END DO
1153	ENDIF
1154	!
1155	ENDIF
1156	!
1157	END SUBROUTINE reflux_sshu
1158
1159	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1160	!!----------------------------------------------------------------------
1161	!! * ROUTINE updatevb2b *
1162	!!----------------------------------------------------------------------
1163	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1164	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1165	LOGICAL , INTENT(in ) :: before
1166	!!
1167	INTEGER :: ji, jj
1168	REAL(wp) :: zrhox, za1, zcor
1169	!!---------------------------------------------
1170	!
1171	IF( before ) THEN
1172	zrhox = Agrif_Rhox()
1173	DO jj=j1,j2
1174	DO ji=i1,i2
1175	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1176	END DO
1177	END DO
1178	tabres = zrhox * tabres
1179	ELSE
1180	!
1181	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1182	!
1183	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1184	DO jj=j1,j2
1185	DO ji=i1,i2
1186	zcor=tabres(ji,jj) - vb2_b(ji,jj)
1187	! Update time integrated fluxes also in case of multiply nested grids:
1188	vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor
1189	! Update corrective fluxes:
1190	vn_bf(ji,jj) = vn_bf(ji,jj) + zcor
1191	! Update half step back fluxes:
1192	vb2_b(ji,jj) = tabres(ji,jj)
1193	END DO
1194	END DO
1195	ENDIF
1196	!
1197	END SUBROUTINE updatevb2b
1198
1199	SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1200	!!---------------------------------------------
1201	!! * ROUTINE reflux_sshv *
1202	!!---------------------------------------------
1203	INTEGER, INTENT(in) :: i1, i2, j1, j2
1204	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1205	LOGICAL, INTENT(in) :: before
1206	INTEGER, INTENT(in) :: nb, ndir
1207	!!
1208	LOGICAL :: southern_side, northern_side
1209	INTEGER :: ji, jj
1210	REAL(wp) :: zrhox, za1, zcor
1211	!!---------------------------------------------
1212	!
1213	IF (before) THEN
1214	zrhox = Agrif_Rhox()
1215	DO jj=j1,j2
1216	DO ji=i1,i2
1217	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1218	END DO
1219	END DO
1220	tabres = zrhox * tabres
1221	ELSE
1222	!
1223	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1224	!
1225	southern_side = (nb == 2).AND.(ndir == 1)
1226	northern_side = (nb == 2).AND.(ndir == 2)
1227	!
1228	IF (southern_side) THEN
1229	DO ji=i1,i2
1230	zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1))
1231	ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor
1232	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor
1233	END DO
1234	ENDIF
1235	IF (northern_side) THEN
1236	DO ji=i1,i2
1237	zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2))
1238	ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
1239	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor
1240	END DO
1241	ENDIF
1242	!
1243	ENDIF
1244	!
1245	END SUBROUTINE reflux_sshv
1246
1247	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1248	!
1249	! ====>>>>>>>>>> currently not used
1250	!
1251	!!----------------------------------------------------------------------
1252	!! * ROUTINE updateT *
1253	!!----------------------------------------------------------------------
1254	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
1255	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
1256	LOGICAL , INTENT(in ) :: before
1257	!!
1258	INTEGER :: ji,jj,jk
1259	REAL(wp) :: ztemp
1260	!!----------------------------------------------------------------------
1261
1262	IF (before) THEN
1263	DO jk=k1,k2
1264	DO jj=j1,j2
1265	DO ji=i1,i2
1266	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
1267	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1268	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1269	END DO
1270	END DO
1271	END DO
1272	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
1273	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1274	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1275	ELSE
1276	DO jk=k1,k2
1277	DO jj=j1,j2
1278	DO ji=i1,i2
1279	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1280	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
1281	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
1282	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
1283	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1284	print *,'CORR = ',ztemp-1.
1285	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
1286	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
1287	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1288	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1289	END IF
1290	END DO
1291	END DO
1292	END DO
1293	ENDIF
1294	!
1295	END SUBROUTINE update_scales
1296
1297
1298	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1299	!!----------------------------------------------------------------------
1300	!! * ROUTINE updateen *
1301	!!----------------------------------------------------------------------
1302	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1303	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1304	LOGICAL , INTENT(in ) :: before
1305	!!----------------------------------------------------------------------
1306	!
1307	IF( before ) THEN
1308	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1309	ELSE
1310	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1311	ENDIF
1312	!
1313	END SUBROUTINE updateEN
1314
1315
1316	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1317	!!----------------------------------------------------------------------
1318	!! * ROUTINE updateavt *
1319	!!----------------------------------------------------------------------
1320	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1321	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1322	LOGICAL , INTENT(in ) :: before
1323	!!----------------------------------------------------------------------
1324	!
1325	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1326	ELSE ; avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1327	ENDIF
1328	!
1329	END SUBROUTINE updateAVT
1330
1331
1332	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1333	!!---------------------------------------------
1334	!! * ROUTINE updateavm *
1335	!!----------------------------------------------------------------------
1336	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1337	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1338	LOGICAL , INTENT(in ) :: before
1339	!!----------------------------------------------------------------------
1340	!
1341	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1342	ELSE ; avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1343	ENDIF
1344	!
1345	END SUBROUTINE updateAVM
1346
1347	#if ! defined key_qco
1348	SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1349	!!---------------------------------------------
1350	!! * ROUTINE updatee3t *
1351	!!---------------------------------------------
1352	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1353	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1354	LOGICAL, INTENT(in) :: before
1355	!
1356	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1357	INTEGER :: ji,jj,jk
1358	REAL(wp) :: zcoef
1359	!!---------------------------------------------
1360	!
1361	IF (.NOT.before) THEN
1362	!
1363	ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))
1364	!
1365	! Update e3t from ssh (z* case only)
1366	DO jk = 1, jpkm1
1367	DO jj=j1,j2
1368	DO ji=i1,i2
1369	ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) &
1370	& *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1371	END DO
1372	END DO
1373	END DO
1374	!
1375	! 1) Updates at BEFORE time step:
1376	! -------------------------------
1377	!
1378	! Save "old" scale factor (prior update) for subsequent asselin correction
1379	! of prognostic variables
1380	e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a)
1381
1382	! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace...
1383	! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a)
1384
1385	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
1386	DO jk = 1, jpkm1
1387	DO jj=j1,j2
1388	DO ji=i1,i2
1389	e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) &
1390	& + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
1391	END DO
1392	END DO
1393	END DO
1394	!
1395	e3w (i1:i2,j1:j2,1,Kbb_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb_a) - e3t_0(i1:i2,j1:j2,1)
1396	gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp
1397	gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a)
1398	!
1399	DO jk = 2, jpk
1400	DO jj = j1,j2
1401	DO ji = i1,i2
1402	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1403	e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &
1404	& ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) &
1405	& + 0.5_wp * tmask(ji,jj,jk) * &
1406	& ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) )
1407	gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
1408	gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) &
1409	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a))
1410	END DO
1411	END DO
1412	END DO
1413	!
1414	ENDIF
1415	!
1416	! 2) Updates at NOW time step:
1417	! ----------------------------
1418	!
1419	! Update vertical scale factor at T-points:
1420	e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1)
1421	!
1422	! Update total depth:
1423	ht(i1:i2,j1:j2) = 0._wp
1424	DO jk = 1, jpkm1
1425	ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk)
1426	END DO
1427	!
1428	! Update vertical scale factor at W-points and depths:
1429	e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)
1430	gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a)
1431	gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp
1432	gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
1433	!
1434	DO jk = 2, jpk
1435	DO jj = j1,j2
1436	DO ji = i1,i2
1437	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1438	e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) &
1439	& + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) )
1440	gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
1441	gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) &
1442	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a))
1443	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
1444	END DO
1445	END DO
1446	END DO
1447	!
1448	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1449	e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a)
1450	e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a)
1451	gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a)
1452	gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
1453	ENDIF
1454	!
1455	DEALLOCATE(ptab)
1456	ENDIF
1457	!
1458	END SUBROUTINE updatee3t
1459	#endif
1460
1461	#else
1462	!!----------------------------------------------------------------------
1463	!! Empty module no AGRIF zoom
1464	!!----------------------------------------------------------------------
1465	CONTAINS
1466	SUBROUTINE agrif_oce_update_empty
1467	WRITE(,) 'agrif_oce_update : You should not have seen this print! error?'
1468	END SUBROUTINE agrif_oce_update_empty
1469	#endif
1470
1471	!!======================================================================
1472	END MODULE agrif_oce_update
1473

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

Download in other formats: