New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

agrif_oce_update.F90 in NEMO/trunk/src/NST – NEMO

Context Navigation

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

Last change on this file since 14064 was 14064, checked in by ayoung, 3 years ago
Merging ticket #2506 into trunk.
Property svn:keywords set to `Id`
File size: 59.3 KB

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

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

Download in other formats: