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/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90 @ 11603

Last change on this file since 11603 was 11603, checked in by jchanut, 5 years ago
#2222, 1) correct time interpolation of barotropic velocities in corners. 2) Clean remapping module and enable remapping several variables at the same time. At this stage, vertical remapping doesn't change VORTEX results with an identical vertical grid ONLY in one way mode and a linearized free surface (within truncature errors).
Property svn:keywords set to `Id`
File size: 54.7 KB

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

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

Download in other formats: