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

Last change on this file since 11590 was 11590, checked in by jchanut, 5 years ago
#2222: 1) create remapping module (vremap) and integration of D. Engwirda piecewise polynomial recontruction package (PPR_LIB cpp key). 2) Various bug corrections with key_vertical activated.
Property svn:keywords set to `Id`
File size: 54.8 KB

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

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

Download in other formats: