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

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

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

Download in other formats: