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

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

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

Context Navigation

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

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

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

Download in other formats: