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 @ 10030

Last change on this file since 10030 was 10030, checked in by gm, 6 years ago
#1911 (ENHANCE-04): RK3 branch - step II.3 remove e3uw_$ e3vw_$, except e3.w_0 and use only after e3 in dyn/trazdf
Property svn:keywords set to `Id`
File size: 57.3 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	hu_a (:,:) = hu_n (:,:)
261	hv_a (:,:) = hv_n (:,:)
262	!
263	! !* NOW fields :
264	CALL ssh2e3_now ! set: ht, hu, hv, r1_hu, r1_hv
265	! ! e3t, e3u, e3v, e3f (from 1 to jpkm1)
266	! ! e3w, gdept_n, gdepw_n, gde3w_n (from 1 to jpk )
267
268	! !* BEFORE fields :
269	IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN
270	!
271	CALL ssh2e3_before ! set: hu, hv, r1_hu, r1_hv
272	! ! e3t, e3u, e3v (from 1 to jpkm1)
273	! ! e3w (from 1 to jpk )
274	ENDIF
275	!
276	END SUBROUTINE dom_vvl_update_UVF
277
278	# if defined key_vertical
279
280	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
281	!!----------------------------------------------------------------------
282	!! * ROUTINE updateT *
283	!!----------------------------------------------------------------------
284	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
285	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
286	LOGICAL, INTENT(in) :: before
287	!!
288	INTEGER :: ji,jj,jk,jn
289	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
290	REAL(wp) :: h_in(k1:k2)
291	REAL(wp) :: h_out(1:jpk)
292	INTEGER :: N_in, N_out
293	REAL(wp) :: zrho_xy, h_diff
294	REAL(wp) :: tabin(k1:k2,n1:n2)
295	!!----------------------------------------------------------------------
296	!
297	IF (before) THEN
298	AGRIF_SpecialValue = -999._wp
299	zrho_xy = Agrif_rhox() * Agrif_rhoy()
300	DO jn = n1, n2-1
301	DO jk = k1, k2
302	DO jj = j1, j2
303	DO ji = i1, i2
304	tabres(ji,jj,jk,jn) = tmask(ji,jj,jk) * (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
305	+ (tmask(ji,jj,jk)-1) * 999._wp
306	END DO
307	END DO
308	END DO
309	END DO
310	DO jk = k1, k2
311	DO jj = j1, j2
312	DO ji = i1, i2
313	tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
314	& + (tmask(ji,jj,jk)-1) * 999._wp
315	END DO
316	END DO
317	END DO
318	ELSE
319	tabres_child(:,:,:,:) = 0.
320	AGRIF_SpecialValue = 0._wp
321	DO jj = j1 , j2
322	DO ji = i1, i2
323	N_in = 0
324	DO jk = k1, k2 !k2 = jpk of child grid
325	IF ( tabres(ji,jj,jk,n2) == 0 ) EXIT
326	N_in = N_in + 1
327	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
328	h_in (N_in) = tabres(ji,jj,jk,n2)
329	END DO
330	N_out = 0
331	DO jk = 1, jpk ! jpk of parent grid
332	IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
333	N_out = N_out + 1
334	h_out(N_out) = e3t_n(ji,jj,jk)
335	END DO
336	IF (N_in > 0) THEN !Remove this?
337	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
338	IF (h_diff < -1.e-4) THEN
339	print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
340	print *,h_in(1:N_in)
341	print *,h_out(1:N_out)
342	STOP
343	ENDIF
344	DO jn = n1, n2-1
345	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)
346	END DO
347	ENDIF
348	END DO
349	END DO
350
351	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
352	DO jn = n1, n2-1
353	DO jk = 1, jpk
354	DO jj = j1, j2
355	DO ji = i1, i2
356	IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN
357	tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &
358	& + rn_atfp * ( tabres_child(ji,jj,jk,jn) - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
359	ENDIF
360	END DO
361	END DO
362	END DO
363	END DO
364	ENDIF
365	DO jn = n1, n2-1
366	DO jk = 1, jpk
367	DO jj = j1, j2
368	DO ji = i1, i2
369	IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN
370	tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
371	END IF
372	END DO
373	END DO
374	END DO
375	END DO
376	ENDIF
377	!
378	END SUBROUTINE updateTS
379
380	# else
381
382	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
383	!!----------------------------------------------------------------------
384	!! * ROUTINE ROUTINE updateT *
385	!!----------------------------------------------------------------------
386	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
387	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
388	LOGICAL, INTENT(in) :: before
389	!
390	INTEGER :: ji,jj,jk,jn
391	REAL(wp) :: ztb, ztnu, ztno
392	!!----------------------------------------------------------------------
393	!
394	IF (before) THEN
395	DO jn = 1,jpts
396	DO jk=k1,k2
397	DO jj=j1,j2
398	DO ji=i1,i2
399	!> jc tmp
400	tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
401	! tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk)
402	!< jc tmp
403	END DO
404	END DO
405	END DO
406	END DO
407	ELSE
408	!> jc tmp
409	DO jn = 1,jpts
410	tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
411	& * tmask(i1:i2,j1:j2,k1:k2)
412	END DO
413	!< jc tmp
414	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
415	DO jn = 1,jpts
416	DO jk = k1, k2
417	DO jj = j1, j2
418	DO ji = i1, i2
419	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
420	ztb = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
421	ztnu = tabres(ji,jj,jk,jn)
422	ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
423	tsb(ji,jj,jk,jn) = ( ztb + rn_atfp * ( ztnu - ztno) ) / e3t_b(ji,jj,jk) * tmask(ji,jj,jk)
424	ENDIF
425	END DO
426	END DO
427	END DO
428	END DO
429	ENDIF
430	DO jn = 1,jpts
431	DO jk=k1,k2
432	DO jj=j1,j2
433	DO ji=i1,i2
434	IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
435	tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
436	END IF
437	END DO
438	END DO
439	END DO
440	END DO
441	!
442	IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
443	tsb(i1:i2,j1:j2,k1:k2,1:jpts) = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
444	ENDIF
445	!
446	ENDIF
447	!
448	END SUBROUTINE updateTS
449
450	# endif
451
452	# if defined key_vertical
453
454	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
455	!!----------------------------------------------------------------------
456	!! * ROUTINE updateu *
457	!!----------------------------------------------------------------------
458	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
459	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
460	LOGICAL , INTENT(in ) :: before
461	!
462	INTEGER :: ji, jj, jk
463	REAL(wp):: zrhoy
464	! VERTICAL REFINEMENT BEGIN
465	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
466	REAL(wp) :: h_in(k1:k2)
467	REAL(wp) :: h_out(1:jpk)
468	INTEGER :: N_in, N_out
469	REAL(wp) :: h_diff, excess, thick
470	REAL(wp) :: tabin(k1:k2)
471	! VERTICAL REFINEMENT END
472	!!----------------------------------------------------------------------
473	!
474	IF( before ) THEN
475	zrhoy = Agrif_Rhoy()
476	AGRIF_SpecialValue = -999._wp
477	DO jk=k1,k2
478	DO jj=j1,j2
479	DO ji=i1,i2
480	tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) &
481	+ (umask(ji,jj,jk)-1)*999._wp
482	tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) &
483	+ (umask(ji,jj,jk)-1)*999._wp
484	END DO
485	END DO
486	END DO
487	ELSE
488	tabres_child(:,:,:) = 0.
489	AGRIF_SpecialValue = 0._wp
490	DO jj=j1,j2
491	DO ji=i1,i2
492	N_in = 0
493	h_in(:) = 0._wp
494	tabin(:) = 0._wp
495	DO jk=k1,k2 !k2=jpk of child grid
496	IF( tabres(ji,jj,jk,2) < -900) EXIT
497	N_in = N_in + 1
498	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
499	h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
500	END DO
501	N_out = 0
502	DO jk=1,jpk
503	IF (umask(ji,jj,jk) == 0) EXIT
504	N_out = N_out + 1
505	h_out(N_out) = e3u_n(ji,jj,jk)
506	END DO
507	IF (N_in * N_out > 0) THEN
508	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
509	IF (h_diff < -1.e-4) THEN
510	!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
511	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
512	excess = 0._wp
513	DO jk=N_in,1,-1
514	thick = MIN(-1*h_diff, h_in(jk))
515	excess = excess + tabin(jk)thicke2u(ji,jj)
516	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
517	h_diff = h_diff + thick
518	IF ( h_diff == 0) THEN
519	N_in = jk
520	h_in(jk) = h_in(jk) - thick
521	EXIT
522	ENDIF
523	END DO
524	ENDIF
525	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)
526	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
527	ENDIF
528	END DO
529	END DO
530
531	DO jk = 1, jpk
532	DO jj = j1, j2
533	DO ji = i1, i2
534	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler) ) THEN ! Add asselin part
535	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
536	ub(ji,jj,jk) = ub(ji,jj,jk) + rn_atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
537	ENDIF
538	un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
539	END DO
540	END DO
541	END DO
542	ENDIF
543	!
544	END SUBROUTINE updateu
545
546	#else
547
548	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
549	!!----------------------------------------------------------------------
550	!! * ROUTINE updateu *
551	!!----------------------------------------------------------------------
552	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
553	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
554	LOGICAL , INTENT(in ) :: before
555	!
556	INTEGER :: ji, jj, jk
557	REAL(wp) :: zrhoy, zub, zunu, zuno
558	!!----------------------------------------------------------------------
559	!
560	IF( before ) THEN
561	zrhoy = Agrif_Rhoy()
562	DO jk = k1, k2
563	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)
564	END DO
565	ELSE
566	DO jk=k1,k2
567	DO jj=j1,j2
568	DO ji=i1,i2
569	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)
570	!
571	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
572	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
573	zub = ub(ji,jj,jk) * e3u_b(ji,jj,jk) ! fse3t_b prior update should be used
574	zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
575	zunu = tabres(ji,jj,jk,1)
576	ub(ji,jj,jk) = ( zub + rn_atfp * ( zunu - zuno) ) / e3u_b(ji,jj,jk) * umask(ji,jj,jk)
577	ENDIF
578	!
579	un(ji,jj,jk) = tabres(ji,jj,jk,1) / e3u_n(ji,jj,jk) * umask(ji,jj,jk)
580	END DO
581	END DO
582	END DO
583	!
584	IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
585	!!gm IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
586	ub(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
587	ENDIF
588	!
589	ENDIF
590	!
591	END SUBROUTINE updateu
592
593	# endif
594
595	SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
596	!!----------------------------------------------------------------------
597	!! * ROUTINE correct_u_bdy *
598	!!----------------------------------------------------------------------
599	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
600	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
601	LOGICAL , INTENT(in ) :: before
602	INTEGER , INTENT(in ) :: nb, ndir
603	!!
604	LOGICAL :: western_side, eastern_side
605	INTEGER :: jj, jk
606	REAL(wp):: zcor
607	!!----------------------------------------------------------------------
608	!
609	IF( .NOT.before ) THEN
610	!
611	western_side = (nb == 1).AND.(ndir == 1)
612	eastern_side = (nb == 1).AND.(ndir == 2)
613	!
614	IF (western_side) THEN
615	DO jj=j1,j2
616	zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
617	un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
618	DO jk=1,jpkm1
619	un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
620	END DO
621	END DO
622	ENDIF
623	!
624	IF (eastern_side) THEN
625	DO jj=j1,j2
626	zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
627	un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
628	DO jk=1,jpkm1
629	un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
630	END DO
631	END DO
632	ENDIF
633	!
634	ENDIF
635	!
636	END SUBROUTINE correct_u_bdy
637
638	# if defined key_vertical
639
640	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
641	!!----------------------------------------------------------------------
642	!! * ROUTINE updatev *
643	!!----------------------------------------------------------------------
644	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
645	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
646	LOGICAL , INTENT(in ) :: before
647	!
648	INTEGER :: ji, jj, jk
649	REAL(wp) :: zrhox
650	! VERTICAL REFINEMENT BEGIN
651	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
652	REAL(wp) :: h_in(k1:k2)
653	REAL(wp) :: h_out(1:jpk)
654	INTEGER :: N_in, N_out
655	REAL(wp) :: h_diff, excess, thick
656	REAL(wp) :: tabin(k1:k2)
657	! VERTICAL REFINEMENT END
658	!!----------------------------------------------------------------------
659	!
660	IF( before ) THEN
661	zrhox = Agrif_Rhox()
662	AGRIF_SpecialValue = -999._wp
663	DO jk=k1,k2
664	DO jj=j1,j2
665	DO ji=i1,i2
666	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
667	+ (vmask(ji,jj,jk)-1)*999._wp
668	tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
669	+ (vmask(ji,jj,jk)-1)*999._wp
670	END DO
671	END DO
672	END DO
673	ELSE
674	tabres_child(:,:,:) = 0.
675	AGRIF_SpecialValue = 0._wp
676	DO jj=j1,j2
677	DO ji=i1,i2
678	N_in = 0
679	DO jk=k1,k2
680	IF (tabres(ji,jj,jk,2) < -900) EXIT
681	N_in = N_in + 1
682	tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
683	h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
684	END DO
685	N_out = 0
686	DO jk=1,jpk
687	IF (vmask(ji,jj,jk) == 0) EXIT
688	N_out = N_out + 1
689	h_out(N_out) = e3v_n(ji,jj,jk)
690	END DO
691	IF (N_in * N_out > 0) THEN
692	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
693	IF (h_diff < -1.e-4) then
694	!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
695	!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
696	excess = 0._wp
697	DO jk=N_in,1,-1
698	thick = MIN(-1*h_diff, h_in(jk))
699	excess = excess + tabin(jk)thicke2u(ji,jj)
700	tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
701	h_diff = h_diff + thick
702	IF ( h_diff == 0) THEN
703	N_in = jk
704	h_in(jk) = h_in(jk) - thick
705	EXIT
706	ENDIF
707	END DO
708	ENDIF
709	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)
710	tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
711	ENDIF
712	END DO
713	END DO
714
715	DO jk=1,jpk
716	DO jj=j1,j2
717	DO ji=i1,i2
718	!
719	IF( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
720	!!gm IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
721	vb(ji,jj,jk) = vb(ji,jj,jk) + rn_atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
722	ENDIF
723	!
724	vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
725	END DO
726	END DO
727	END DO
728	ENDIF
729	!
730	END SUBROUTINE updatev
731
732	# else
733
734	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
735	!!----------------------------------------------------------------------
736	!! * ROUTINE updatev *
737	!!----------------------------------------------------------------------
738	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
739	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
740	LOGICAL , INTENT(in ) :: before
741	!
742	INTEGER :: ji, jj, jk
743	REAL(wp) :: zrhox, zvb, zvnu, zvno
744	!!----------------------------------------------------------------------
745	!
746	IF (before) THEN
747	zrhox = Agrif_Rhox()
748	DO jk=k1,k2
749	DO jj=j1,j2
750	DO ji=i1,i2
751	tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
752	END DO
753	END DO
754	END DO
755	ELSE
756	DO jk=k1,k2
757	DO jj=j1,j2
758	DO ji=i1,i2
759	tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
760	!
761	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
762	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
763	zvb = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
764	zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
765	zvnu = tabres(ji,jj,jk,1)
766	vb(ji,jj,jk) = ( zvb + rn_atfp * ( zvnu - zvno) ) / e3v_b(ji,jj,jk) * vmask(ji,jj,jk)
767	ENDIF
768	!
769	vn(ji,jj,jk) = tabres(ji,jj,jk,1) / e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
770	END DO
771	END DO
772	END DO
773	!
774	IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
775	!!gm IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
776	vb(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
777	ENDIF
778	!
779	ENDIF
780	!
781	END SUBROUTINE updatev
782
783	# endif
784
785	SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
786	!!----------------------------------------------------------------------
787	!! * ROUTINE correct_v_bdy *
788	!!----------------------------------------------------------------------
789	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
790	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
791	LOGICAL , INTENT(in ) :: before
792	INTEGER , INTENT(in) :: nb, ndir
793	!!
794	LOGICAL :: southern_side, northern_side
795	!
796	INTEGER :: ji, jk
797	REAL(wp) :: zcor
798	!!----------------------------------------------------------------------
799	!
800	IF( .NOT.before ) THEN
801	!
802	southern_side = (nb == 2).AND.(ndir == 1)
803	northern_side = (nb == 2).AND.(ndir == 2)
804	!
805	IF (southern_side) THEN
806	DO ji=i1,i2
807	zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
808	vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
809	DO jk=1,jpkm1
810	vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
811	END DO
812	END DO
813	ENDIF
814	!
815	IF (northern_side) THEN
816	DO ji=i1,i2
817	zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
818	vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
819	DO jk=1,jpkm1
820	vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
821	END DO
822	END DO
823	ENDIF
824	!
825	ENDIF
826	!
827	END SUBROUTINE correct_v_bdy
828
829
830	SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
831	!!----------------------------------------------------------------------
832	!! * ROUTINE updateu2d *
833	!!----------------------------------------------------------------------
834	INTEGER , INTENT(in ) :: i1, i2, j1, j2
835	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
836	LOGICAL , INTENT(in ) :: before
837	!
838	INTEGER :: ji, jj, jk
839	REAL(wp) :: zrhoy
840	REAL(wp) :: zcorr
841	!!----------------------------------------------------------------------
842	!
843	IF( before ) THEN
844	zrhoy = Agrif_Rhoy()
845	DO jj=j1,j2
846	DO ji=i1,i2
847	tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
848	END DO
849	END DO
850	ELSE
851	DO jj=j1,j2
852	DO ji=i1,i2
853	tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj)
854	!
855	! Update "now" 3d velocities:
856	spgu(ji,jj) = 0._wp
857	DO jk=1,jpkm1
858	spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
859	END DO
860	!
861	zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
862	DO jk=1,jpkm1
863	un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)
864	END DO
865	!
866	! Update barotropic velocities:
867	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
868	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
869	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
870	zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
871	ub_b(ji,jj) = ub_b(ji,jj) + rn_atfp * zcorr * umask(ji,jj,1)
872	END IF
873	ENDIF
874	un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
875	!
876	! Correct "before" velocities to hold correct bt component:
877	spgu(ji,jj) = 0.e0
878	DO jk=1,jpkm1
879	spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
880	END DO
881	!
882	zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
883	DO jk=1,jpkm1
884	ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)
885	END DO
886	!
887	END DO
888	END DO
889	!
890	IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
891	!!gm IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
892	ub_b(i1:i2,j1:j2) = un_b(i1:i2,j1:j2)
893	ENDIF
894	ENDIF
895	!
896	END SUBROUTINE updateu2d
897
898
899	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
900	!!----------------------------------------------------------------------
901	!! * ROUTINE updatev2d *
902	!!----------------------------------------------------------------------
903	INTEGER , INTENT(in ) :: i1, i2, j1, j2
904	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
905	LOGICAL , INTENT(in ) :: before
906	!
907	INTEGER :: ji, jj, jk
908	REAL(wp) :: zrhox, zcorr
909	!!----------------------------------------------------------------------
910	!
911	IF( before ) THEN
912	zrhox = Agrif_Rhox()
913	DO jj=j1,j2
914	DO ji=i1,i2
915	tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj)
916	END DO
917	END DO
918	ELSE
919	DO jj=j1,j2
920	DO ji=i1,i2
921	tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj)
922	!
923	! Update "now" 3d velocities:
924	spgv(ji,jj) = 0.e0
925	DO jk=1,jpkm1
926	spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
927	END DO
928	!
929	zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
930	DO jk=1,jpkm1
931	vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)
932	END DO
933	!
934	! Update barotropic velocities:
935	IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN
936	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part
937	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
938	zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
939	vb_b(ji,jj) = vb_b(ji,jj) + rn_atfp * zcorr * vmask(ji,jj,1)
940	END IF
941	ENDIF
942	vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
943	!
944	! Correct "before" velocities to hold correct bt component:
945	spgv(ji,jj) = 0.e0
946	DO jk=1,jpkm1
947	spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
948	END DO
949	!
950	zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
951	DO jk=1,jpkm1
952	vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)
953	END DO
954	!
955	END DO
956	END DO
957	!
958	IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
959	!!gm IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
960	vb_b(i1:i2,j1:j2) = vn_b(i1:i2,j1:j2)
961	ENDIF
962	!
963	ENDIF
964	!
965	END SUBROUTINE updatev2d
966
967
968	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
969	!!----------------------------------------------------------------------
970	!! * ROUTINE updateSSH *
971	!!----------------------------------------------------------------------
972	INTEGER , INTENT(in ) :: i1, i2, j1, j2
973	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
974	LOGICAL , INTENT(in ) :: before
975	!
976	INTEGER :: ji, jj
977	!!----------------------------------------------------------------------
978	!
979	IF( before ) THEN
980	DO jj = j1, j2
981	DO ji = i1, i2
982	tabres(ji,jj) = ssh(ji,jj,Nnn)
983	END DO
984	END DO
985	ELSE
986	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN
987	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
988	DO jj = j1, j2
989	DO ji = i1, i2
990	ssh(ji,jj,Nbb) = ssh(ji,jj,Nbb) + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Nnn) ) * tmask(ji,jj,1)
991	END DO
992	END DO
993	ENDIF
994	!
995	DO jj = j1, j2
996	DO ji = i1, i2
997	ssh(ji,jj,Nnn) = tabres(ji,jj) * tmask(ji,jj,1)
998	END DO
999	END DO
1000	!
1001	IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN
1002	!!gm IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1003	ssh(i1:i2,j1:j2,Nbb) = ssh(i1:i2,j1:j2,Nnn)
1004	ENDIF
1005	!
1006	ENDIF
1007	!
1008	END SUBROUTINE updateSSH
1009
1010
1011	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1012	!!----------------------------------------------------------------------
1013	!! * ROUTINE updateub2b *
1014	!!----------------------------------------------------------------------
1015	INTEGER , INTENT(in) :: i1, i2, j1, j2
1016	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1017	LOGICAL , INTENT(in) :: before
1018	!!
1019	INTEGER :: ji, jj
1020	REAL(wp) :: zrhoy, za1, zcor
1021	!!---------------------------------------------
1022	!
1023	IF (before) THEN
1024	zrhoy = Agrif_Rhoy()
1025	DO jj=j1,j2
1026	DO ji=i1,i2
1027	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1028	END DO
1029	END DO
1030	tabres = zrhoy * tabres
1031	ELSE
1032	!
1033	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1034	!
1035	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1036	DO jj=j1,j2
1037	DO ji=i1,i2
1038	zcor=tabres(ji,jj) - ub2_b(ji,jj)
1039	! Update time integrated fluxes also in case of multiply nested grids:
1040	ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor
1041	! Update corrective fluxes:
1042	un_bf(ji,jj) = un_bf(ji,jj) + zcor
1043	! Update half step back fluxes:
1044	ub2_b(ji,jj) = tabres(ji,jj)
1045	END DO
1046	END DO
1047	ENDIF
1048	!
1049	END SUBROUTINE updateub2b
1050
1051
1052	SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1053	!!----------------------------------------------------------------------
1054	!! * ROUTINE reflux_sshu *
1055	!!----------------------------------------------------------------------
1056	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1057	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1058	LOGICAL , INTENT(in ) :: before
1059	INTEGER , INTENT(in ) :: nb, ndir
1060	!
1061	LOGICAL :: western_side, eastern_side
1062	INTEGER :: ji, jj
1063	REAL(wp):: zrhoy, za1, zcor
1064	!!----------------------------------------------------------------------
1065	!
1066	IF (before) THEN
1067	zrhoy = Agrif_Rhoy()
1068	DO jj=j1,j2
1069	DO ji=i1,i2
1070	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1071	END DO
1072	END DO
1073	tabres = zrhoy * tabres
1074	ELSE
1075	!
1076	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1077	!
1078	western_side = (nb == 1).AND.(ndir == 1)
1079	eastern_side = (nb == 1).AND.(ndir == 2)
1080	!
1081	IF ( western_side ) THEN
1082	DO jj=j1,j2
1083	zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
1084	ssh(i1 ,jj,Nnn) = ssh(i1 ,jj,Nnn) + zcor
1085	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) ssh(i1 ,jj,Nbb) = ssh(i1 ,jj,Nbb) + rn_atfp * zcor
1086	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1 ,jj) = sshb(i1 ,jj) + rn_atfp * zcor
1087	END DO
1088	ENDIF
1089	IF ( eastern_side ) THEN
1090	DO jj=j1,j2
1091	zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1092	ssh(i2+1,jj,Nnn) = ssh(i2+1,jj,Nnn) + zcor
1093	IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) ssh(i2+1,jj,Nbb) = ssh(i2+1,jj,Nbb) + rn_atfp * zcor
1094	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor
1095	END DO
1096	ENDIF
1097	!
1098	ENDIF
1099	!
1100	END SUBROUTINE reflux_sshu
1101
1102
1103	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1104	!!----------------------------------------------------------------------
1105	!! * ROUTINE updatevb2b *
1106	!!----------------------------------------------------------------------
1107	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1108	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1109	LOGICAL , INTENT(in ) :: before
1110	!
1111	INTEGER :: ji, jj
1112	REAL(wp) :: zrhox, za1, zcor
1113	!!---------------------------------------------------------------------
1114	!
1115	IF( before ) THEN
1116	zrhox = Agrif_Rhox()
1117	DO jj=j1,j2
1118	DO ji=i1,i2
1119	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1120	END DO
1121	END DO
1122	tabres = zrhox * tabres
1123	ELSE
1124	!
1125	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1126	!
1127	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1128	DO jj=j1,j2
1129	DO ji=i1,i2
1130	zcor=tabres(ji,jj) - vb2_b(ji,jj)
1131	! Update time integrated fluxes also in case of multiply nested grids:
1132	vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor
1133	! Update corrective fluxes:
1134	vn_bf(ji,jj) = vn_bf(ji,jj) + zcor
1135	! Update half step back fluxes:
1136	vb2_b(ji,jj) = tabres(ji,jj)
1137	END DO
1138	END DO
1139	ENDIF
1140	!
1141	END SUBROUTINE updatevb2b
1142
1143
1144	SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1145	!!----------------------------------------------------------------------
1146	!! * ROUTINE reflux_sshv *
1147	!!----------------------------------------------------------------------
1148	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1149	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1150	LOGICAL , INTENT(in ) :: before
1151	INTEGER , INTENT(in ) :: nb, ndir
1152	!!
1153	LOGICAL :: southern_side, northern_side
1154	INTEGER :: ji, jj
1155	REAL(wp) :: zrhox, za1, zcor
1156	!!----------------------------------------------------------------------
1157	!
1158	IF (before) THEN
1159	zrhox = Agrif_Rhox()
1160	DO jj=j1,j2
1161	DO ji=i1,i2
1162	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
1163	END DO
1164	END DO
1165	tabres = zrhox * tabres
1166	ELSE
1167	!
1168	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1169	!
1170	southern_side = (nb == 2).AND.(ndir == 1)
1171	northern_side = (nb == 2).AND.(ndir == 2)
1172	!
1173	IF (southern_side) THEN
1174	DO ji=i1,i2
1175	zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * ( vb2_b(ji,j1)-tabres(ji,j1) )
1176	ssh(ji,j1 ,Nnn) = ssh(ji,j1 ,Nnn) + zcor
1177	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) ssh(ji,j1 ,Nbb) = ssh(ji,j1,Nbb) + rn_atfp * zcor
1178	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1 ) = sshb(ji,j1) + rn_atfp * zcor
1179	END DO
1180	ENDIF
1181	IF (northern_side) THEN
1182	DO ji=i1,i2
1183	zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * ( vb2_b(ji,j2)-tabres(ji,j2) )
1184	ssh(ji,j2+1,Nnn) = ssh(ji,j2+1,Nnn) + zcor
1185	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) ssh(ji,j2+1,Nbb) = ssh(ji,j2+1,Nbb) + rn_atfp * zcor
1186	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor
1187	END DO
1188	ENDIF
1189	!
1190	ENDIF
1191	!
1192	END SUBROUTINE reflux_sshv
1193
1194
1195	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1196	!!----------------------------------------------------------------------
1197	!! * ROUTINE updateT *
1198	!
1199	! ====>>>>>>>>>> currently not used
1200	!
1201	!!----------------------------------------------------------------------
1202	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
1203	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
1204	LOGICAL , INTENT(in ) :: before
1205	!!
1206	INTEGER :: ji,jj,jk
1207	REAL(wp) :: ztemp
1208	!!----------------------------------------------------------------------
1209
1210	IF (before) THEN
1211	DO jk=k1,k2
1212	DO jj=j1,j2
1213	DO ji=i1,i2
1214	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
1215	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1216	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1217	END DO
1218	END DO
1219	END DO
1220	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
1221	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1222	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1223	ELSE
1224	DO jk=k1,k2
1225	DO jj=j1,j2
1226	DO ji=i1,i2
1227	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1228	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
1229	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
1230	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
1231	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1232	print *,'CORR = ',ztemp-1.
1233	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
1234	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
1235	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1236	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1237	END IF
1238	END DO
1239	END DO
1240	END DO
1241	ENDIF
1242	!
1243	END SUBROUTINE update_scales
1244
1245
1246	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1247	!!----------------------------------------------------------------------
1248	!! * ROUTINE updateen *
1249	!!----------------------------------------------------------------------
1250	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1251	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1252	LOGICAL , INTENT(in ) :: before
1253	!!----------------------------------------------------------------------
1254	!
1255	IF( before ) THEN
1256	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1257	ELSE
1258	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1259	ENDIF
1260	!
1261	END SUBROUTINE updateEN
1262
1263
1264	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1265	!!----------------------------------------------------------------------
1266	!! * ROUTINE updateavt *
1267	!!----------------------------------------------------------------------
1268	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1269	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1270	LOGICAL , INTENT(in ) :: before
1271	!!----------------------------------------------------------------------
1272	!
1273	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1274	ELSE ; avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1275	ENDIF
1276	!
1277	END SUBROUTINE updateAVT
1278
1279
1280	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1281	!!----------------------------------------------------------------------
1282	!! * ROUTINE updateavm *
1283	!!----------------------------------------------------------------------
1284	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1285	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1286	LOGICAL , INTENT(in ) :: before
1287	!!----------------------------------------------------------------------
1288	!
1289	IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1290	ELSE ; avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1291	ENDIF
1292	!
1293	END SUBROUTINE updateAVM
1294
1295
1296	SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1297	!!----------------------------------------------------------------------
1298	!! * ROUTINE updatee3t *
1299	!!----------------------------------------------------------------------
1300	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1301	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1302	LOGICAL, INTENT(in) :: before
1303	!
1304	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1305	INTEGER :: ji,jj,jk
1306	REAL(wp) :: zcoef
1307	!!---------------------------------------------
1308	!
1309	IF (.NOT.before) THEN
1310	!
1311	ALLOCATE( ptab(i1:i2,j1:j2,1:jpk) )
1312	!
1313	! Update e3t from ssh (z* case only)
1314	DO jk = 1, jpkm1
1315	DO jj = j1, j2
1316	DO ji = i1, i2
1317	ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Nnn) * r1_ht_0(ji,jj) *tmask(ji,jj,jk) )
1318	END DO
1319	END DO
1320	END DO
1321	!
1322	! 1) Updates at BEFORE time step:
1323	! -------------------------------
1324	!
1325	! Save "old" scale factor (prior update) for subsequent asselin correction
1326	! of prognostic variables
1327	e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1328
1329	! One should also save e3t_b, but lacking of workspace...
1330	! hdivn(i1:i2,j1:j2,1:jpkm1) = e3t_b(i1:i2,j1:j2,1:jpkm1)
1331
1332	IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN
1333	!!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1334	DO jk = 1, jpkm1
1335	DO jj=j1,j2
1336	DO ji=i1,i2
1337	e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) + rn_atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1338	END DO
1339	END DO
1340	END DO
1341	!
1342	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)
1343	gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1344	gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1345	!
1346	DO jk = 2, jpk
1347	DO jj = j1,j2
1348	DO ji = i1,i2
1349	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1350	e3w_b(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &
1351	& ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) &
1352	& + 0.5_wp * tmask(ji,jj,jk) * &
1353	& ( e3t_b(ji,jj,jk ) - e3t_0(ji,jj,jk ) )
1354	gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1355	gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &
1356	& + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk))
1357	END DO
1358	END DO
1359	END DO
1360	!
1361	ENDIF
1362	!
1363	! 2) Updates at NOW time step:
1364	! ----------------------------
1365	!
1366	! Update vertical scale factor at T-points:
1367	e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1368	!
1369	! Update total depth:
1370	ht_n(i1:i2,j1:j2) = 0._wp
1371	DO jk = 1, jpkm1
1372	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)
1373	END DO
1374	!
1375	! Update vertical scale factor at W-points and depths:
1376	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)
1377	gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1378	gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1379	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
1380	!
1381	DO jk = 2, jpk
1382	DO jj = j1,j2
1383	DO ji = i1,i2
1384	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1385	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) ) &
1386	& + 0.5_wp * tmask(ji,jj,jk) * ( e3t_n(ji,jj,jk ) - e3t_0(ji,jj,jk ) )
1387	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1388	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &
1389	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
1390	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
1391	END DO
1392	END DO
1393	END DO
1394	!
1395	IF ( l_1st_euler .AND. Agrif_Nb_Step()==0 ) THEN
1396	!!gm IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1397	e3t_b (i1:i2,j1:j2,1:jpk) = e3t_n (i1:i2,j1:j2,1:jpk)
1398	e3w_b (i1:i2,j1:j2,1:jpk) = e3w_n (i1:i2,j1:j2,1:jpk)
1399	gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1400	gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1401	ENDIF
1402	!
1403	DEALLOCATE(ptab)
1404	ENDIF
1405	!
1406	END SUBROUTINE updatee3t
1407
1408	#else
1409	!!----------------------------------------------------------------------
1410	!! Empty module no AGRIF zoom
1411	!!----------------------------------------------------------------------
1412	CONTAINS
1413	SUBROUTINE agrif_oce_update_empty
1414	WRITE(,) 'agrif_oce_update : You should not have seen this print! error?'
1415	END SUBROUTINE agrif_oce_update_empty
1416	#endif
1417
1418	!!======================================================================
1419	END MODULE agrif_oce_update
1420

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

Download in other formats: