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_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_update.F90 @ 14048

Last change on this file since 14048 was 13937, checked in by jchanut, 4 years ago

#2222, added:

Divergence conserving interpolation for transport (Balsara's scheme)
Possibility to shift barotropic interface inside the domain (currently disabled - needed to match volume within sponge)
PPR library as default vertical interpolation scheme
Trilinear interpolation over partial bottom cells

Property svn:keywords set to Id

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

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

Download in other formats: