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/trunk/src/NST – NEMO

Context Navigation

source: NEMO/trunk/src/NST/agrif_oce_update.F90 @ 14053

Last change on this file since 14053 was 14053, checked in by techene, 3 years ago
#2385 added to the trunk
Property svn:keywords set to `Id`
File size: 59.1 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	CALL Agrif_Update_qco
222	CALL Agrif_ParentGrid_To_ChildGrid()
223	#endif
224	!
225	END SUBROUTINE Agrif_Update_vvl
226
227
228	#if defined key_qco
229	SUBROUTINE Agrif_Update_qco
230	!!---------------------------------------------
231	!! * ROUTINE dom_Update_qco *
232	!!---------------------------------------------
233	!
234	! Save arrays prior update (needed for asselin correction)
235	r3t(:,:,Krhs_a) = r3t(:,:,Kmm_a)
236	r3u(:,:,Krhs_a) = r3u(:,:,Kmm_a)
237	r3v(:,:,Krhs_a) = r3v(:,:,Kmm_a)
238
239	! Update r3x arrays from updated ssh
240	CALL dom_qco_zgr( Kbb_a, Kmm_a )
241	!
242	END SUBROUTINE Agrif_Update_qco
243	#endif
244
245
246	#if ! defined key_qco
247	SUBROUTINE dom_vvl_update_UVF
248	!!---------------------------------------------
249	!! * ROUTINE dom_vvl_update_UVF *
250	!!---------------------------------------------
251	!!
252	INTEGER :: jk
253	REAL(wp):: zcoef
254	!!---------------------------------------------
255	IF (lwp.AND.lk_agrif_debug) Write(,) 'Finalize e3 on grid Number', &
256	& Agrif_Fixed(), 'Step', Agrif_Nb_Step()
257
258	! Save "old" scale factor (prior update) for subsequent asselin correction
259	! of prognostic variables
260	! -----------------------
261	!
262	e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a)
263	e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a)
264	! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a)
265	! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a)
266	hu(:,:,Krhs_a) = hu(:,:,Kmm_a)
267	hv(:,:,Krhs_a) = hv(:,:,Kmm_a)
268
269	! 1) NOW fields
270	!--------------
271
272	! Vertical scale factor interpolations
273	! ------------------------------------
274	CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) , 'U' )
275	CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) , 'V' )
276	CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) , 'F' )
277
278	CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' )
279	CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' )
280
281	! Update total depths:
282	! --------------------
283	hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points
284	hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points
285	DO jk = 1, jpkm1
286	hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)
287	hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk)
288	END DO
289	! ! Inverse of the local depth
290	r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) )
291	r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) )
292
293
294	! 2) BEFORE fields:
295	!------------------
296	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
297	!
298	! Vertical scale factor interpolations
299	! ------------------------------------
300	CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a), 'U' )
301	CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a), 'V' )
302
303	CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' )
304	CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' )
305
306	! Update total depths:
307	! --------------------
308	hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points
309	hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points
310	DO jk = 1, jpkm1
311	hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk)
312	hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk)
313	END DO
314	! ! Inverse of the local depth
315	r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) )
316	r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) )
317	ENDIF
318	!
319	END SUBROUTINE dom_vvl_update_UVF
320	#endif
321
322	#if defined key_vertical
323
324	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
325	!!----------------------------------------------------------------------
326	!! * ROUTINE updateT *
327	!!---------------------------------------------
328	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
329	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
330	LOGICAL, INTENT(in) :: before
331	!!
332	INTEGER :: ji,jj,jk,jn
333	INTEGER :: N_in, N_out
334	REAL(wp) :: ztb, ztnu, ztno
335	REAL(wp) :: h_in(k1:k2)
336	REAL(wp) :: h_out(1:jpk)
337	REAL(wp) :: tabin(k1:k2,1:jpts)
338	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child
339	!!---------------------------------------------
340	!
341	IF (before) THEN
342	!jc_alt
343	! AGRIF_SpecialValue = -999._wp
344	DO jn = n1,n2-1
345	DO jk=k1,k2
346	DO jj=j1,j2
347	DO ji=i1,i2
348	!jc_alt
349	! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
350	! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp
351	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
352	END DO
353	END DO
354	END DO
355	END DO
356	DO jk=k1,k2
357	DO jj=j1,j2
358	DO ji=i1,i2
359	!jc_alt
360	! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
361	! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp
362	tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
363	END DO
364	END DO
365	END DO
366	ELSE
367	tabres_child(:,:,:,:) = 0._wp
368	AGRIF_SpecialValue = 0._wp
369	DO jj=j1,j2
370	DO ji=i1,i2
371	N_in = 0
372	DO jk=k1,k2 !k2 = jpk of child grid
373	! jc_alt
374	! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT
375	IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT
376	N_in = N_in + 1
377	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
378	h_in(N_in) = tabres(ji,jj,jk,n2)
379	ENDDO
380	N_out = 0
381	DO jk=1,jpk ! jpk of parent grid
382	IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
383	N_out = N_out + 1
384	h_out(N_out) = e3t(ji,jj,jk,Kmm_a)
385	ENDDO
386	IF (N_in*N_out > 0) THEN !Remove this?
387	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)
388	ENDIF
389	ENDDO
390	ENDDO
391
392	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
393	! Add asselin part
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	ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
400	ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
401	ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
402	ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &
403	& * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
404	ENDIF
405	END DO
406	END DO
407	END DO
408	END DO
409	ENDIF
410	DO jn = 1,jpts
411	DO jk = 1, jpkm1
412	DO jj = j1, j2
413	DO ji = i1, i2
414	IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
415	ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
416	END IF
417	END DO
418	END DO
419	END DO
420	END DO
421	!
422	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
423	ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a)
424	ENDIF
425	ENDIF
426	!
427	END SUBROUTINE updateTS
428
429	# else
430
431	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
432	!!---------------------------------------------
433	!! * ROUTINE updateT *
434	!!---------------------------------------------
435	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
436	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
437	LOGICAL, INTENT(in) :: before
438	!!
439	INTEGER :: ji,jj,jk,jn
440	REAL(wp) :: ztb, ztnu, ztno
441	!!---------------------------------------------
442	!
443	IF (before) THEN
444	DO jn = 1,jpts
445	DO jk=k1,k2
446	DO jj=j1,j2
447	DO ji=i1,i2
448	!> jc tmp
449	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
450	! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
451	!< jc tmp
452	END DO
453	END DO
454	END DO
455	END DO
456	ELSE
457	!> jc tmp
458	DO jn = 1,jpts
459	tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
460	& * tmask(i1:i2,j1:j2,k1:k2)
461	ENDDO
462	!< jc tmp
463	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
464	! Add asselin part
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	ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
471	ztnu = tabres(ji,jj,jk,jn)
472	ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
473	ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &
474	& * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
475	ENDIF
476	END DO
477	END DO
478	END DO
479	END DO
480	ENDIF
481	DO jn = 1,jpts
482	DO jk=k1,k2
483	DO jj=j1,j2
484	DO ji=i1,i2
485	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
486	ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
487	END IF
488	END DO
489	END DO
490	END DO
491	END DO
492	!
493	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
494	ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)
495	ENDIF
496	!
497	ENDIF
498	!
499	END SUBROUTINE updateTS
500
501	# endif
502
503	# if defined key_vertical
504
505	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
506	!!---------------------------------------------
507	!! * ROUTINE updateu *
508	!!---------------------------------------------
509	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
510	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
511	LOGICAL , INTENT(in ) :: before
512	!
513	INTEGER :: ji, jj, jk
514	REAL(wp):: zrhoy, zub, zunu, zuno
515	! VERTICAL REFINEMENT BEGIN
516	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
517	REAL(wp) :: h_in(k1:k2)
518	REAL(wp) :: h_out(1:jpk)
519	INTEGER :: N_in, N_out
520	REAL(wp) :: h_diff, excess, thick
521	REAL(wp) :: tabin(k1:k2)
522	! VERTICAL REFINEMENT END
523	!!---------------------------------------------
524	!
525	IF( before ) THEN
526	zrhoy = Agrif_Rhoy()
527	!jc_alt
528	! AGRIF_SpecialValue = -999._wp
529	DO jk=k1,k2
530	DO jj=j1,j2
531	DO ji=i1,i2
532	!jc_alt
533	! 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) &
534	! & + (umask(ji,jj,jk)-1._wp)*999._wp
535	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)
536	!jc_alt
537	! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) &
538	! & + (umask(ji,jj,jk)-1._wp)*999._wp
539	tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
540	END DO
541	END DO
542	END DO
543	ELSE
544	tabres_child(:,:,:) = 0.
545	AGRIF_SpecialValue = 0._wp
546	DO jj=j1,j2
547	DO ji=i1,i2
548	N_in = 0
549	h_in(:) = 0._wp
550	tabin(:) = 0._wp
551	DO jk=k1,k2 !k2=jpk of child grid
552	!jc_alt
553	! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT
554	IF( tabres(ji,jj,jk,2) == 0.) EXIT
555	N_in = N_in + 1
556	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
557	h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
558	ENDDO
559	N_out = 0
560	DO jk=1,jpk
561	IF (umask(ji,jj,jk) == 0) EXIT
562	N_out = N_out + 1
563	h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
564	ENDDO
565	IF (N_in * N_out > 0) THEN
566	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
567	excess = 0._wp
568	IF (h_diff < -1.e-4) THEN
569	!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
570	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
571	DO jk=N_in,1,-1
572	thick = MIN(-1*h_diff, h_in(jk))
573	excess = excess + tabin(jk)thicke2u(ji,jj)
574	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
575	h_diff = h_diff + thick
576	IF ( h_diff == 0) THEN
577	N_in = jk
578	h_in(jk) = h_in(jk) - thick
579	EXIT
580	ENDIF
581	ENDDO
582	ENDIF
583	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)
584	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
585	ENDIF
586	ENDDO
587	ENDDO
588	!
589	DO jk=1,jpk
590	DO jj=j1,j2
591	DO ji=i1,i2
592	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
593	zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
594	zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
595	zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a)
596	uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &
597	& * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
598	ENDIF
599	!
600	uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
601	END DO
602	END DO
603	END DO
604	!
605	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
606	uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a)
607	ENDIF
608	!
609	ENDIF
610	!
611	END SUBROUTINE updateu
612
613	#else
614
615	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
616	!!---------------------------------------------
617	!! * ROUTINE updateu *
618	!!---------------------------------------------
619	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
620	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
621	LOGICAL , INTENT(in ) :: before
622	!
623	INTEGER :: ji, jj, jk
624	REAL(wp) :: zrhoy, zub, zunu, zuno
625	!!---------------------------------------------
626	!
627	IF( before ) THEN
628	zrhoy = Agrif_Rhoy()
629	DO jk = k1, k2
630	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)
631	END DO
632	ELSE
633	DO jk=k1,k2
634	DO jj=j1,j2
635	DO ji=i1,i2
636	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)
637	!
638	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
639	zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
640	zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
641	zunu = tabres(ji,jj,jk,1)
642	uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &
643	& * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
644	ENDIF
645	!
646	uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a)
647	END DO
648	END DO
649	END DO
650	!
651	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
652	uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a)
653	ENDIF
654	!
655	ENDIF
656	!
657	END SUBROUTINE updateu
658
659	# endif
660
661	SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
662	!!---------------------------------------------
663	!! * ROUTINE correct_u_bdy *
664	!!---------------------------------------------
665	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
666	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
667	LOGICAL , INTENT(in ) :: before
668	INTEGER , INTENT(in) :: nb, ndir
669	!!
670	LOGICAL :: western_side, eastern_side
671	!
672	INTEGER :: jj, jk
673	REAL(wp) :: zcor
674	!!---------------------------------------------
675	!
676	IF( .NOT.before ) THEN
677	!
678	western_side = (nb == 1).AND.(ndir == 1)
679	eastern_side = (nb == 1).AND.(ndir == 2)
680	!
681	IF (western_side) THEN
682	DO jj=j1,j2
683	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)
684	uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor
685	DO jk=1,jpkm1
686	uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk)
687	END DO
688	END DO
689	ENDIF
690	!
691	IF (eastern_side) THEN
692	DO jj=j1,j2
693	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)
694	uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor
695	DO jk=1,jpkm1
696	uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk)
697	END DO
698	END DO
699	ENDIF
700	!
701	ENDIF
702	!
703	END SUBROUTINE correct_u_bdy
704
705	# if defined key_vertical
706
707	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
708	!!---------------------------------------------
709	!! * ROUTINE updatev *
710	!!---------------------------------------------
711	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
712	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
713	LOGICAL , INTENT(in ) :: before
714	!
715	INTEGER :: ji, jj, jk
716	REAL(wp) :: zrhox, zvb, zvnu, zvno
717	! VERTICAL REFINEMENT BEGIN
718	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
719	REAL(wp) :: h_in(k1:k2)
720	REAL(wp) :: h_out(1:jpk)
721	INTEGER :: N_in, N_out
722	REAL(wp) :: h_diff, excess, thick
723	REAL(wp) :: tabin(k1:k2)
724	! VERTICAL REFINEMENT END
725	!!---------------------------------------------
726	!
727	IF( before ) THEN
728	zrhox = Agrif_Rhox()
729	!jc_alt
730	! AGRIF_SpecialValue = -999._wp
731	DO jk=k1,k2
732	DO jj=j1,j2
733	DO ji=i1,i2
734	!jc_alt
735	! 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) &
736	! & + (vmask(ji,jj,jk)-1._wp) * 999._wp
737	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)
738	!jc_alt
739	! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) &
740	! & + (vmask(ji,jj,jk)-1._wp) * 999._wp
741	tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk)
742	END DO
743	END DO
744	END DO
745	ELSE
746	tabres_child(:,:,:) = 0.
747	AGRIF_SpecialValue = 0._wp
748	DO jj=j1,j2
749	DO ji=i1,i2
750	N_in = 0
751	DO jk=k1,k2
752	!jc_alt
753	! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT
754	IF (tabres(ji,jj,jk,2) == 0) EXIT
755	N_in = N_in + 1
756	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
757	h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
758	ENDDO
759	N_out = 0
760	DO jk=1,jpk
761	IF (vmask(ji,jj,jk) == 0) EXIT
762	N_out = N_out + 1
763	h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
764	ENDDO
765	IF (N_in * N_out > 0) THEN
766	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
767	excess = 0._wp
768	IF (h_diff < -1.e-4) then
769	!Even if bathy at T points match it's possible for the V points to be deeper in the child grid.
770	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
771	DO jk=N_in,1,-1
772	thick = MIN(-1*h_diff, h_in(jk))
773	excess = excess + tabin(jk)thicke2u(ji,jj)
774	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
775	h_diff = h_diff + thick
776	IF ( h_diff == 0) THEN
777	N_in = jk
778	h_in(jk) = h_in(jk) - thick
779	EXIT
780	ENDIF
781	ENDDO
782	ENDIF
783	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)
784	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
785	ENDIF
786	ENDDO
787	ENDDO
788	!
789	DO jk=1,jpkm1
790	DO jj=j1,j2
791	DO ji=i1,i2
792	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
793	zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
794	zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
795	zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a)
796	vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &
797	& * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
798	ENDIF
799	!
800	vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
801	END DO
802	END DO
803	END DO
804	!
805	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
806	vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a)
807	ENDIF
808	!
809	ENDIF
810	!
811	END SUBROUTINE updatev
812
813	# else
814
815	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
816	!!---------------------------------------------
817	!! * ROUTINE updatev *
818	!!---------------------------------------------
819	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
820	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
821	LOGICAL , INTENT(in ) :: before
822	!
823	INTEGER :: ji, jj, jk
824	REAL(wp) :: zrhox, zvb, zvnu, zvno
825	!!---------------------------------------------
826	!
827	IF (before) THEN
828	zrhox = Agrif_Rhox()
829	DO jk=k1,k2
830	DO jj=j1,j2
831	DO ji=i1,i2
832	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
833	END DO
834	END DO
835	END DO
836	ELSE
837	DO jk=k1,k2
838	DO jj=j1,j2
839	DO ji=i1,i2
840	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
841	!
842	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
843	zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
844	zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
845	zvnu = tabres(ji,jj,jk,1)
846	vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &
847	& * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
848	ENDIF
849	!
850	vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a)
851	END DO
852	END DO
853	END DO
854	!
855	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
856	vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a)
857	ENDIF
858	!
859	ENDIF
860	!
861	END SUBROUTINE updatev
862
863	# endif
864
865	SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
866	!!---------------------------------------------
867	!! * ROUTINE correct_v_bdy *
868	!!---------------------------------------------
869	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
870	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
871	LOGICAL , INTENT(in ) :: before
872	INTEGER , INTENT(in) :: nb, ndir
873	!!
874	LOGICAL :: southern_side, northern_side
875	!
876	INTEGER :: ji, jk
877	REAL(wp) :: zcor
878	!!---------------------------------------------
879	!
880	IF( .NOT.before ) THEN
881	!
882	southern_side = (nb == 2).AND.(ndir == 1)
883	northern_side = (nb == 2).AND.(ndir == 2)
884	!
885	IF (southern_side) THEN
886	DO ji=i1,i2
887	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)
888	vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor
889	DO jk=1,jpkm1
890	vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk)
891	END DO
892	END DO
893	ENDIF
894	!
895	IF (northern_side) THEN
896	DO ji=i1,i2
897	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)
898	vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor
899	DO jk=1,jpkm1
900	vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk)
901	END DO
902	END DO
903	ENDIF
904	!
905	ENDIF
906	!
907	END SUBROUTINE correct_v_bdy
908
909
910	SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
911	!!----------------------------------------------------------------------
912	!! * ROUTINE updateu2d *
913	!!----------------------------------------------------------------------
914	INTEGER , INTENT(in ) :: i1, i2, j1, j2
915	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
916	LOGICAL , INTENT(in ) :: before
917	!!
918	INTEGER :: ji, jj, jk
919	REAL(wp) :: zrhoy
920	REAL(wp) :: zcorr
921	!!---------------------------------------------
922	!
923	IF( before ) THEN
924	zrhoy = Agrif_Rhoy()
925	DO jj=j1,j2
926	DO ji=i1,i2
927	tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj)
928	END DO
929	END DO
930	ELSE
931	DO jj=j1,j2
932	DO ji=i1,i2
933	tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj)
934	!
935	! Update "now" 3d velocities:
936	spgu(ji,jj) = 0._wp
937	DO jk=1,jpkm1
938	spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
939	END DO
940	!
941	zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)
942	DO jk=1,jpkm1
943	uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)
944	END DO
945	!
946	! Update barotropic velocities:
947	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
948	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
949	zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a)
950	uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1)
951	END IF
952	ENDIF
953	uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1)
954	!
955	! Correct "before" velocities to hold correct bt component:
956	spgu(ji,jj) = 0.e0
957	DO jk=1,jpkm1
958	spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
959	END DO
960	!
961	zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a)
962	DO jk=1,jpkm1
963	uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)
964	END DO
965	!
966	END DO
967	END DO
968	!
969	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
970	uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a)
971	ENDIF
972	ENDIF
973	!
974	END SUBROUTINE updateu2d
975
976
977	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
978	!!----------------------------------------------------------------------
979	!! * ROUTINE updatev2d *
980	!!----------------------------------------------------------------------
981	INTEGER , INTENT(in ) :: i1, i2, j1, j2
982	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
983	LOGICAL , INTENT(in ) :: before
984	!
985	INTEGER :: ji, jj, jk
986	REAL(wp) :: zrhox, zcorr
987	!!----------------------------------------------------------------------
988	!
989	IF( before ) THEN
990	zrhox = Agrif_Rhox()
991	DO jj=j1,j2
992	DO ji=i1,i2
993	tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj)
994	END DO
995	END DO
996	ELSE
997	DO jj=j1,j2
998	DO ji=i1,i2
999	tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj)
1000	!
1001	! Update "now" 3d velocities:
1002	spgv(ji,jj) = 0.e0
1003	DO jk=1,jpkm1
1004	spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
1005	END DO
1006	!
1007	zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)
1008	DO jk=1,jpkm1
1009	vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)
1010	END DO
1011	!
1012	! Update barotropic velocities:
1013	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
1014	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
1015	zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a)
1016	vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1)
1017	END IF
1018	ENDIF
1019	vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1)
1020	!
1021	! Correct "before" velocities to hold correct bt component:
1022	spgv(ji,jj) = 0.e0
1023	DO jk=1,jpkm1
1024	spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
1025	END DO
1026	!
1027	zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a)
1028	DO jk=1,jpkm1
1029	vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)
1030	END DO
1031	!
1032	END DO
1033	END DO
1034	!
1035	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1036	vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a)
1037	ENDIF
1038	!
1039	ENDIF
1040	!
1041	END SUBROUTINE updatev2d
1042
1043
1044	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
1045	!!----------------------------------------------------------------------
1046	!! * ROUTINE updateSSH *
1047	!!----------------------------------------------------------------------
1048	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1049	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1050	LOGICAL , INTENT(in ) :: before
1051	!!
1052	INTEGER :: ji, jj
1053	!!----------------------------------------------------------------------
1054	!
1055	IF( before ) THEN
1056	DO jj=j1,j2
1057	DO ji=i1,i2
1058	tabres(ji,jj) = ssh(ji,jj,Kmm_a)
1059	END DO
1060	END DO
1061	ELSE
1062	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
1063	DO jj=j1,j2
1064	DO ji=i1,i2
1065	ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) &
1066	& + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
1067	END DO
1068	END DO
1069	ENDIF
1070	!
1071	DO jj=j1,j2
1072	DO ji=i1,i2
1073	ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1)
1074	END DO
1075	END DO
1076	!
1077	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1078	ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a)
1079	ENDIF
1080	!
1081
1082	ENDIF
1083	!
1084	END SUBROUTINE updateSSH
1085
1086
1087	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1088	!!----------------------------------------------------------------------
1089	!! * ROUTINE updateub2b *
1090	!!----------------------------------------------------------------------
1091	INTEGER , INTENT(in) :: i1, i2, j1, j2
1092	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1093	LOGICAL , INTENT(in) :: before
1094	!!
1095	INTEGER :: ji, jj
1096	REAL(wp) :: zrhoy, za1, zcor
1097	!!---------------------------------------------
1098	!
1099	IF (before) THEN
1100	zrhoy = Agrif_Rhoy()
1101	DO jj=j1,j2
1102	DO ji=i1,i2
1103	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1104	END DO
1105	END DO
1106	tabres = zrhoy * tabres
1107	ELSE
1108	!
1109	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1110	!
1111	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1112	DO jj=j1,j2
1113	DO ji=i1,i2
1114	zcor=tabres(ji,jj) - ub2_b(ji,jj)
1115	! Update time integrated fluxes also in case of multiply nested grids:
1116	ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor
1117	! Update corrective fluxes:
1118	un_bf(ji,jj) = un_bf(ji,jj) + zcor
1119	! Update half step back fluxes:
1120	ub2_b(ji,jj) = tabres(ji,jj)
1121	END DO
1122	END DO
1123	ENDIF
1124	!
1125	END SUBROUTINE updateub2b
1126
1127	SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1128	!!---------------------------------------------
1129	!! * ROUTINE reflux_sshu *
1130	!!---------------------------------------------
1131	INTEGER, INTENT(in) :: i1, i2, j1, j2
1132	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1133	LOGICAL, INTENT(in) :: before
1134	INTEGER, INTENT(in) :: nb, ndir
1135	!!
1136	LOGICAL :: western_side, eastern_side
1137	INTEGER :: ji, jj
1138	REAL(wp) :: zrhoy, za1, zcor
1139	!!---------------------------------------------
1140	!
1141	IF (before) THEN
1142	zrhoy = Agrif_Rhoy()
1143	DO jj=j1,j2
1144	DO ji=i1,i2
1145	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1146	END DO
1147	END DO
1148	tabres = zrhoy * tabres
1149	ELSE
1150	!
1151	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1152	!
1153	western_side = (nb == 1).AND.(ndir == 1)
1154	eastern_side = (nb == 1).AND.(ndir == 2)
1155	!
1156	IF (western_side) THEN
1157	DO jj=j1,j2
1158	zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
1159	ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor
1160	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + rn_atfp * zcor
1161	END DO
1162	ENDIF
1163	IF (eastern_side) THEN
1164	DO jj=j1,j2
1165	zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1166	ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
1167	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor
1168	END DO
1169	ENDIF
1170	!
1171	ENDIF
1172	!
1173	END SUBROUTINE reflux_sshu
1174
1175	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1176	!!----------------------------------------------------------------------
1177	!! * ROUTINE updatevb2b *
1178	!!----------------------------------------------------------------------
1179	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1180	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1181	LOGICAL , INTENT(in ) :: before
1182	!!
1183	INTEGER :: ji, jj
1184	REAL(wp) :: zrhox, za1, zcor
1185	!!---------------------------------------------
1186	!
1187	IF( before ) THEN
1188	zrhox = Agrif_Rhox()
1189	DO jj=j1,j2
1190	DO ji=i1,i2
1191	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1192	END DO
1193	END DO
1194	tabres = zrhox * tabres
1195	ELSE
1196	!
1197	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1198	!
1199	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1200	DO jj=j1,j2
1201	DO ji=i1,i2
1202	zcor=tabres(ji,jj) - vb2_b(ji,jj)
1203	! Update time integrated fluxes also in case of multiply nested grids:
1204	vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor
1205	! Update corrective fluxes:
1206	vn_bf(ji,jj) = vn_bf(ji,jj) + zcor
1207	! Update half step back fluxes:
1208	vb2_b(ji,jj) = tabres(ji,jj)
1209	END DO
1210	END DO
1211	ENDIF
1212	!
1213	END SUBROUTINE updatevb2b
1214
1215	SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1216	!!---------------------------------------------
1217	!! * ROUTINE reflux_sshv *
1218	!!---------------------------------------------
1219	INTEGER, INTENT(in) :: i1, i2, j1, j2
1220	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1221	LOGICAL, INTENT(in) :: before
1222	INTEGER, INTENT(in) :: nb, ndir
1223	!!
1224	LOGICAL :: southern_side, northern_side
1225	INTEGER :: ji, jj
1226	REAL(wp) :: zrhox, za1, zcor
1227	!!---------------------------------------------
1228	!
1229	IF (before) THEN
1230	zrhox = Agrif_Rhox()
1231	DO jj=j1,j2
1232	DO ji=i1,i2
1233	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1234	END DO
1235	END DO
1236	tabres = zrhox * tabres
1237	ELSE
1238	!
1239	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1240	!
1241	southern_side = (nb == 2).AND.(ndir == 1)
1242	northern_side = (nb == 2).AND.(ndir == 2)
1243	!
1244	IF (southern_side) THEN
1245	DO ji=i1,i2
1246	zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1))
1247	ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor
1248	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor
1249	END DO
1250	ENDIF
1251	IF (northern_side) THEN
1252	DO ji=i1,i2
1253	zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2))
1254	ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
1255	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor
1256	END DO
1257	ENDIF
1258	!
1259	ENDIF
1260	!
1261	END SUBROUTINE reflux_sshv
1262
1263	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1264	!
1265	! ====>>>>>>>>>> currently not used
1266	!
1267	!!----------------------------------------------------------------------
1268	!! * ROUTINE updateT *
1269	!!----------------------------------------------------------------------
1270	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
1271	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
1272	LOGICAL , INTENT(in ) :: before
1273	!!
1274	INTEGER :: ji,jj,jk
1275	REAL(wp) :: ztemp
1276	!!----------------------------------------------------------------------
1277
1278	IF (before) THEN
1279	DO jk=k1,k2
1280	DO jj=j1,j2
1281	DO ji=i1,i2
1282	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
1283	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1284	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1285	END DO
1286	END DO
1287	END DO
1288	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
1289	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1290	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1291	ELSE
1292	DO jk=k1,k2
1293	DO jj=j1,j2
1294	DO ji=i1,i2
1295	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1296	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
1297	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
1298	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
1299	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1300	print *,'CORR = ',ztemp-1.
1301	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
1302	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
1303	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1304	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1305	END IF
1306	END DO
1307	END DO
1308	END DO
1309	ENDIF
1310	!
1311	END SUBROUTINE update_scales
1312
1313
1314	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1315	!!----------------------------------------------------------------------
1316	!! * ROUTINE updateen *
1317	!!----------------------------------------------------------------------
1318	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1319	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1320	LOGICAL , INTENT(in ) :: before
1321	!!----------------------------------------------------------------------
1322	!
1323	IF( before ) THEN
1324	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1325	ELSE
1326	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1327	ENDIF
1328	!
1329	END SUBROUTINE updateEN
1330
1331
1332	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1333	!!----------------------------------------------------------------------
1334	!! * ROUTINE updateavt *
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) = avt_k(i1:i2,j1:j2,k1:k2)
1342	ELSE ; avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1343	ENDIF
1344	!
1345	END SUBROUTINE updateAVT
1346
1347
1348	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1349	!!---------------------------------------------
1350	!! * ROUTINE updateavm *
1351	!!----------------------------------------------------------------------
1352	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1353	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1354	LOGICAL , INTENT(in ) :: before
1355	!!----------------------------------------------------------------------
1356	!
1357	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1358	ELSE ; avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1359	ENDIF
1360	!
1361	END SUBROUTINE updateAVM
1362
1363	#if ! defined key_qco
1364	SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1365	!!---------------------------------------------
1366	!! * ROUTINE updatee3t *
1367	!!---------------------------------------------
1368	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1369	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1370	LOGICAL, INTENT(in) :: before
1371	!
1372	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1373	INTEGER :: ji,jj,jk
1374	REAL(wp) :: zcoef
1375	!!---------------------------------------------
1376	!
1377	IF (.NOT.before) THEN
1378	!
1379	ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))
1380	!
1381	! Update e3t from ssh (z* case only)
1382	DO jk = 1, jpkm1
1383	DO jj=j1,j2
1384	DO ji=i1,i2
1385	ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) &
1386	& *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1387	END DO
1388	END DO
1389	END DO
1390	!
1391	! 1) Updates at BEFORE time step:
1392	! -------------------------------
1393	!
1394	! Save "old" scale factor (prior update) for subsequent asselin correction
1395	! of prognostic variables
1396	e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a)
1397
1398	! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace...
1399	! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a)
1400
1401	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
1402	DO jk = 1, jpkm1
1403	DO jj=j1,j2
1404	DO ji=i1,i2
1405	e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) &
1406	& + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
1407	END DO
1408	END DO
1409	END DO
1410	!
1411	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)
1412	gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp
1413	gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a)
1414	!
1415	DO jk = 2, jpk
1416	DO jj = j1,j2
1417	DO ji = i1,i2
1418	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1419	e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &
1420	& ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) &
1421	& + 0.5_wp * tmask(ji,jj,jk) * &
1422	& ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) )
1423	gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
1424	gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) &
1425	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a))
1426	END DO
1427	END DO
1428	END DO
1429	!
1430	ENDIF
1431	!
1432	! 2) Updates at NOW time step:
1433	! ----------------------------
1434	!
1435	! Update vertical scale factor at T-points:
1436	e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1)
1437	!
1438	! Update total depth:
1439	ht(i1:i2,j1:j2) = 0._wp
1440	DO jk = 1, jpkm1
1441	ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk)
1442	END DO
1443	!
1444	! Update vertical scale factor at W-points and depths:
1445	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)
1446	gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a)
1447	gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp
1448	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
1449	!
1450	DO jk = 2, jpk
1451	DO jj = j1,j2
1452	DO ji = i1,i2
1453	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1454	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) ) &
1455	& + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) )
1456	gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
1457	gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) &
1458	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a))
1459	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
1460	END DO
1461	END DO
1462	END DO
1463	!
1464	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1465	e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a)
1466	e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a)
1467	gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a)
1468	gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
1469	ENDIF
1470	!
1471	DEALLOCATE(ptab)
1472	ENDIF
1473	!
1474	END SUBROUTINE updatee3t
1475	#endif
1476
1477	#else
1478	!!----------------------------------------------------------------------
1479	!! Empty module no AGRIF zoom
1480	!!----------------------------------------------------------------------
1481	CONTAINS
1482	SUBROUTINE agrif_oce_update_empty
1483	WRITE(,) 'agrif_oce_update : You should not have seen this print! error?'
1484	END SUBROUTINE agrif_oce_update_empty
1485	#endif
1486
1487	!!======================================================================
1488	END MODULE agrif_oce_update
1489

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

Download in other formats: