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

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

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

Download in other formats: