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_opa_update.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 9031

Last change on this file since 9031 was 9031, checked in by timgraham, 7 years ago
Resolved AGRIF conflicts
Property svn:keywords set to `Id`
File size: 55.0 KB

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

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

Download in other formats: