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/2018/dev_r9838_ENHANCE04_RK3/src/NST – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_oce_update.F90 @ 10010

Last change on this file since 10010 was 10010, checked in by gm, 6 years ago
#1911 (ENHANCE-04): RK3 branch - step II.1 AGRIF compilation OK for time-level dimension on ssh
Property svn:keywords set to `Id`
File size: 57.5 KB

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

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

Download in other formats: