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

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

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

Download in other formats: