Context Navigation

agrif_oce_update.F90 @ 15574

Last change on this file since 15574 was 15574, checked in by techene, 3 years ago
#2605 #2715 trunk merged into dev_r14318_RK3_stage1
Property svn:keywords set to `Id`
File size: 54.5 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/NST/agrif_oce_update.F90 @ 15574

Download in other formats: