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 @ 13678

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

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

Download in other formats: