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_r12558_HPC-08_epico_Extra_Halo/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_update.F90 @ 13229

Last change on this file since 13229 was 13229, checked in by francesca, 4 years ago
dev_r12558_HPC-08_epico_Extra_Halo: merge with trunk@13218, see #2366
Property svn:keywords set to `Id`
File size: 58.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
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.
88
89	use_sign_north = .TRUE.
90	sign_north = -1.
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.
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.
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	INTEGER :: ji, jj, jk
890	REAL(wp) :: zrhoy
891	REAL(wp) :: zcorr
892	!!---------------------------------------------
893	!
894	IF( before ) THEN
895	zrhoy = Agrif_Rhoy()
896	DO jj=j1,j2
897	DO ji=i1,i2
898	tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj)
899	END DO
900	END DO
901	ELSE
902	DO jj=j1,j2
903	DO ji=i1,i2
904	tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj)
905	!
906	! Update "now" 3d velocities:
907	spgu(ji,jj) = 0._wp
908	DO jk=1,jpkm1
909	spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
910	END DO
911	!
912	zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)
913	DO jk=1,jpkm1
914	uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)
915	END DO
916	!
917	! Update barotropic velocities:
918	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
919	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
920	zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a)
921	uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1)
922	END IF
923	ENDIF
924	uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1)
925	!
926	! Correct "before" velocities to hold correct bt component:
927	spgu(ji,jj) = 0.e0
928	DO jk=1,jpkm1
929	spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
930	END DO
931	!
932	zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a)
933	DO jk=1,jpkm1
934	uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)
935	END DO
936	!
937	END DO
938	END DO
939	!
940	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
941	uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a)
942	ENDIF
943	ENDIF
944	!
945	END SUBROUTINE updateu2d
946
947
948	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
949	!!----------------------------------------------------------------------
950	!! * ROUTINE updatev2d *
951	!!----------------------------------------------------------------------
952	INTEGER , INTENT(in ) :: i1, i2, j1, j2
953	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
954	LOGICAL , INTENT(in ) :: before
955	!
956	INTEGER :: ji, jj, jk
957	REAL(wp) :: zrhox, zcorr
958	!!----------------------------------------------------------------------
959	!
960	IF( before ) THEN
961	zrhox = Agrif_Rhox()
962	DO jj=j1,j2
963	DO ji=i1,i2
964	tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj)
965	END DO
966	END DO
967	ELSE
968	DO jj=j1,j2
969	DO ji=i1,i2
970	tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj)
971	!
972	! Update "now" 3d velocities:
973	spgv(ji,jj) = 0.e0
974	DO jk=1,jpkm1
975	spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
976	END DO
977	!
978	zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)
979	DO jk=1,jpkm1
980	vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)
981	END DO
982	!
983	! Update barotropic velocities:
984	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
985	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
986	zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a)
987	vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1)
988	END IF
989	ENDIF
990	vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1)
991	!
992	! Correct "before" velocities to hold correct bt component:
993	spgv(ji,jj) = 0.e0
994	DO jk=1,jpkm1
995	spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
996	END DO
997	!
998	zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a)
999	DO jk=1,jpkm1
1000	vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)
1001	END DO
1002	!
1003	END DO
1004	END DO
1005	!
1006	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1007	vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a)
1008	ENDIF
1009	!
1010	ENDIF
1011	!
1012	END SUBROUTINE updatev2d
1013
1014
1015	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
1016	!!----------------------------------------------------------------------
1017	!! * ROUTINE updateSSH *
1018	!!----------------------------------------------------------------------
1019	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1020	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1021	LOGICAL , INTENT(in ) :: before
1022	!!
1023	INTEGER :: ji, jj
1024	!!----------------------------------------------------------------------
1025	!
1026	IF( before ) THEN
1027	DO jj=j1,j2
1028	DO ji=i1,i2
1029	tabres(ji,jj) = ssh(ji,jj,Kmm_a)
1030	END DO
1031	END DO
1032	ELSE
1033	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
1034	DO jj=j1,j2
1035	DO ji=i1,i2
1036	ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) &
1037	& + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
1038	END DO
1039	END DO
1040	ENDIF
1041	!
1042	DO jj=j1,j2
1043	DO ji=i1,i2
1044	ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1)
1045	END DO
1046	END DO
1047	!
1048	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1049	ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a)
1050	ENDIF
1051	!
1052
1053	ENDIF
1054	!
1055	END SUBROUTINE updateSSH
1056
1057
1058	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1059	!!----------------------------------------------------------------------
1060	!! * ROUTINE updateub2b *
1061	!!----------------------------------------------------------------------
1062	INTEGER , INTENT(in) :: i1, i2, j1, j2
1063	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1064	LOGICAL , INTENT(in) :: before
1065	!!
1066	INTEGER :: ji, jj
1067	REAL(wp) :: zrhoy, za1, zcor
1068	!!---------------------------------------------
1069	!
1070	IF (before) THEN
1071	zrhoy = Agrif_Rhoy()
1072	DO jj=j1,j2
1073	DO ji=i1,i2
1074	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1075	END DO
1076	END DO
1077	tabres = zrhoy * tabres
1078	ELSE
1079	!
1080	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1081	!
1082	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1083	DO jj=j1,j2
1084	DO ji=i1,i2
1085	zcor=tabres(ji,jj) - ub2_b(ji,jj)
1086	! Update time integrated fluxes also in case of multiply nested grids:
1087	ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor
1088	! Update corrective fluxes:
1089	un_bf(ji,jj) = un_bf(ji,jj) + zcor
1090	! Update half step back fluxes:
1091	ub2_b(ji,jj) = tabres(ji,jj)
1092	END DO
1093	END DO
1094	ENDIF
1095	!
1096	END SUBROUTINE updateub2b
1097
1098	SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1099	!!---------------------------------------------
1100	!! * ROUTINE reflux_sshu *
1101	!!---------------------------------------------
1102	INTEGER, INTENT(in) :: i1, i2, j1, j2
1103	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1104	LOGICAL, INTENT(in) :: before
1105	INTEGER, INTENT(in) :: nb, ndir
1106	!!
1107	LOGICAL :: western_side, eastern_side
1108	INTEGER :: ji, jj
1109	REAL(wp) :: zrhoy, za1, zcor
1110	!!---------------------------------------------
1111	!
1112	IF (before) THEN
1113	zrhoy = Agrif_Rhoy()
1114	DO jj=j1,j2
1115	DO ji=i1,i2
1116	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1117	END DO
1118	END DO
1119	tabres = zrhoy * tabres
1120	ELSE
1121	!
1122	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1123	!
1124	western_side = (nb == 1).AND.(ndir == 1)
1125	eastern_side = (nb == 1).AND.(ndir == 2)
1126	!
1127	IF (western_side) THEN
1128	DO jj=j1,j2
1129	zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
1130	ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor
1131	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + rn_atfp * zcor
1132	END DO
1133	ENDIF
1134	IF (eastern_side) THEN
1135	DO jj=j1,j2
1136	zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1137	ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
1138	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor
1139	END DO
1140	ENDIF
1141	!
1142	ENDIF
1143	!
1144	END SUBROUTINE reflux_sshu
1145
1146	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1147	!!----------------------------------------------------------------------
1148	!! * ROUTINE updatevb2b *
1149	!!----------------------------------------------------------------------
1150	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1151	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1152	LOGICAL , INTENT(in ) :: before
1153	!!
1154	INTEGER :: ji, jj
1155	REAL(wp) :: zrhox, za1, zcor
1156	!!---------------------------------------------
1157	!
1158	IF( before ) THEN
1159	zrhox = Agrif_Rhox()
1160	DO jj=j1,j2
1161	DO ji=i1,i2
1162	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1163	END DO
1164	END DO
1165	tabres = zrhox * tabres
1166	ELSE
1167	!
1168	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1169	!
1170	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1171	DO jj=j1,j2
1172	DO ji=i1,i2
1173	zcor=tabres(ji,jj) - vb2_b(ji,jj)
1174	! Update time integrated fluxes also in case of multiply nested grids:
1175	vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor
1176	! Update corrective fluxes:
1177	vn_bf(ji,jj) = vn_bf(ji,jj) + zcor
1178	! Update half step back fluxes:
1179	vb2_b(ji,jj) = tabres(ji,jj)
1180	END DO
1181	END DO
1182	ENDIF
1183	!
1184	END SUBROUTINE updatevb2b
1185
1186	SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1187	!!---------------------------------------------
1188	!! * ROUTINE reflux_sshv *
1189	!!---------------------------------------------
1190	INTEGER, INTENT(in) :: i1, i2, j1, j2
1191	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1192	LOGICAL, INTENT(in) :: before
1193	INTEGER, INTENT(in) :: nb, ndir
1194	!!
1195	LOGICAL :: southern_side, northern_side
1196	INTEGER :: ji, jj
1197	REAL(wp) :: zrhox, za1, zcor
1198	!!---------------------------------------------
1199	!
1200	IF (before) THEN
1201	zrhox = Agrif_Rhox()
1202	DO jj=j1,j2
1203	DO ji=i1,i2
1204	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1205	END DO
1206	END DO
1207	tabres = zrhox * tabres
1208	ELSE
1209	!
1210	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1211	!
1212	southern_side = (nb == 2).AND.(ndir == 1)
1213	northern_side = (nb == 2).AND.(ndir == 2)
1214	!
1215	IF (southern_side) THEN
1216	DO ji=i1,i2
1217	zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1))
1218	ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor
1219	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor
1220	END DO
1221	ENDIF
1222	IF (northern_side) THEN
1223	DO ji=i1,i2
1224	zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2))
1225	ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
1226	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor
1227	END DO
1228	ENDIF
1229	!
1230	ENDIF
1231	!
1232	END SUBROUTINE reflux_sshv
1233
1234	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1235	!
1236	! ====>>>>>>>>>> currently not used
1237	!
1238	!!----------------------------------------------------------------------
1239	!! * ROUTINE updateT *
1240	!!----------------------------------------------------------------------
1241	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
1242	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
1243	LOGICAL , INTENT(in ) :: before
1244	!!
1245	INTEGER :: ji,jj,jk
1246	REAL(wp) :: ztemp
1247	!!----------------------------------------------------------------------
1248
1249	IF (before) THEN
1250	DO jk=k1,k2
1251	DO jj=j1,j2
1252	DO ji=i1,i2
1253	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
1254	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1255	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1256	END DO
1257	END DO
1258	END DO
1259	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
1260	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1261	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1262	ELSE
1263	DO jk=k1,k2
1264	DO jj=j1,j2
1265	DO ji=i1,i2
1266	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1267	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
1268	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
1269	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
1270	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1271	print *,'CORR = ',ztemp-1.
1272	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
1273	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
1274	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1275	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1276	END IF
1277	END DO
1278	END DO
1279	END DO
1280	ENDIF
1281	!
1282	END SUBROUTINE update_scales
1283
1284
1285	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1286	!!----------------------------------------------------------------------
1287	!! * ROUTINE updateen *
1288	!!----------------------------------------------------------------------
1289	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1290	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1291	LOGICAL , INTENT(in ) :: before
1292	!!----------------------------------------------------------------------
1293	!
1294	IF( before ) THEN
1295	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1296	ELSE
1297	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1298	ENDIF
1299	!
1300	END SUBROUTINE updateEN
1301
1302
1303	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1304	!!----------------------------------------------------------------------
1305	!! * ROUTINE updateavt *
1306	!!----------------------------------------------------------------------
1307	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1308	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1309	LOGICAL , INTENT(in ) :: before
1310	!!----------------------------------------------------------------------
1311	!
1312	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1313	ELSE ; avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1314	ENDIF
1315	!
1316	END SUBROUTINE updateAVT
1317
1318
1319	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1320	!!---------------------------------------------
1321	!! * ROUTINE updateavm *
1322	!!----------------------------------------------------------------------
1323	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1324	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1325	LOGICAL , INTENT(in ) :: before
1326	!!----------------------------------------------------------------------
1327	!
1328	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1329	ELSE ; avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1330	ENDIF
1331	!
1332	END SUBROUTINE updateAVM
1333
1334	SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1335	!!---------------------------------------------
1336	!! * ROUTINE updatee3t *
1337	!!---------------------------------------------
1338	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1339	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1340	LOGICAL, INTENT(in) :: before
1341	!
1342	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1343	INTEGER :: ji,jj,jk
1344	REAL(wp) :: zcoef
1345	!!---------------------------------------------
1346	!
1347	IF (.NOT.before) THEN
1348	!
1349	ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))
1350	!
1351	! Update e3t from ssh (z* case only)
1352	DO jk = 1, jpkm1
1353	DO jj=j1,j2
1354	DO ji=i1,i2
1355	ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) &
1356	& *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1357	END DO
1358	END DO
1359	END DO
1360	!
1361	! 1) Updates at BEFORE time step:
1362	! -------------------------------
1363	!
1364	! Save "old" scale factor (prior update) for subsequent asselin correction
1365	! of prognostic variables
1366	e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a)
1367
1368	! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace...
1369	! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a)
1370
1371	IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN
1372	DO jk = 1, jpkm1
1373	DO jj=j1,j2
1374	DO ji=i1,i2
1375	e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) &
1376	& + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
1377	END DO
1378	END DO
1379	END DO
1380	!
1381	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)
1382	gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp
1383	gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a)
1384	!
1385	DO jk = 2, jpk
1386	DO jj = j1,j2
1387	DO ji = i1,i2
1388	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1389	e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &
1390	& ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) &
1391	& + 0.5_wp * tmask(ji,jj,jk) * &
1392	& ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) )
1393	gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
1394	gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) &
1395	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a))
1396	END DO
1397	END DO
1398	END DO
1399	!
1400	ENDIF
1401	!
1402	! 2) Updates at NOW time step:
1403	! ----------------------------
1404	!
1405	! Update vertical scale factor at T-points:
1406	e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1)
1407	!
1408	! Update total depth:
1409	ht(i1:i2,j1:j2) = 0._wp
1410	DO jk = 1, jpkm1
1411	ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk)
1412	END DO
1413	!
1414	! Update vertical scale factor at W-points and depths:
1415	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)
1416	gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a)
1417	gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp
1418	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
1419	!
1420	DO jk = 2, jpk
1421	DO jj = j1,j2
1422	DO ji = i1,i2
1423	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1424	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) ) &
1425	& + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) )
1426	gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
1427	gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) &
1428	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a))
1429	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
1430	END DO
1431	END DO
1432	END DO
1433	!
1434	IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
1435	e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a)
1436	e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a)
1437	gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a)
1438	gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
1439	ENDIF
1440	!
1441	DEALLOCATE(ptab)
1442	ENDIF
1443	!
1444	END SUBROUTINE updatee3t
1445
1446	#else
1447	!!----------------------------------------------------------------------
1448	!! Empty module no AGRIF zoom
1449	!!----------------------------------------------------------------------
1450	CONTAINS
1451	SUBROUTINE agrif_oce_update_empty
1452	WRITE(,) 'agrif_oce_update : You should not have seen this print! error?'
1453	END SUBROUTINE agrif_oce_update_empty
1454	#endif
1455
1456	!!======================================================================
1457	END MODULE agrif_oce_update
1458

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

Download in other formats: