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

Last change on this file since 11741 was 11741, checked in by jchanut, 4 years ago
#2222: correct definition of parent vertical grid on the child domain to perform vertical interpolation at boundaries. Use additionnal parent depths and number of levels arrays interpolated on the child grid domain to do so. Correction of vertical interpolation of viscosity remains to be done as well as duplication of changes for passive tracers.
Property svn:keywords set to `Id`
File size: 56.5 KB

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

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

Download in other formats: