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

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

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

Download in other formats: