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/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_update.F90 @ 13337

Last change on this file since 13337 was 13337, checked in by jchanut, 4 years ago
#2222, start suppressing key_vertical (add ln_vremap namelist flag)
Property svn:keywords set to `Id`
File size: 53.3 KB

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

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

Download in other formats: