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/2021/dev_r14318_RK3_stage1/src/NST – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/NST/agrif_oce_update.F90 @ 14800

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

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

Download in other formats: