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_r13333_KERNEL-08_techene_gm_HPG_SPG/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/NST/agrif_oce_update.F90 @ 13887

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

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

Download in other formats: