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_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 8901

Last change on this file since 8901 was 8901, checked in by jchanut, 6 years ago
AGRIF: Remove update frequency parameter from namelist - #1965
Property svn:keywords set to `Id`
File size: 41.2 KB

Line
1	#define TWO_WAY /* TWO WAY NESTING */
2	#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
3	#undef VOL_REFLUX /* VOLUME REFLUXING*/
4
5	MODULE agrif_opa_update
6	#if defined key_agrif
7	USE par_oce
8	USE oce
9	USE dom_oce
10	USE agrif_oce
11	USE in_out_manager ! I/O manager
12	USE lib_mpp
13	USE wrk_nemo
14	USE zdf_oce ! vertical physics: ocean variables
15	USE domvvl ! Need interpolation routines
16
17	IMPLICIT NONE
18	PRIVATE
19
20	PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl, Agrif_Update_ssh
21
22	# if defined key_zdftke
23	PUBLIC Agrif_Update_Tke
24	# endif
25	!!----------------------------------------------------------------------
26	!! NEMO/NST 3.6 , NEMO Consortium (2010)
27	!! $Id$
28	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
29	!!----------------------------------------------------------------------
30	CONTAINS
31
32	SUBROUTINE Agrif_Update_Tra( )
33	!!---------------------------------------------
34	!! * ROUTINE Agrif_Update_Tra *
35	!!---------------------------------------------
36	!
37	IF (Agrif_Root()) RETURN
38	!
39	#if defined TWO_WAY
40	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update tracers from grid Number',Agrif_Fixed()
41
42	Agrif_UseSpecialValueInUpdate = .TRUE.
43	Agrif_SpecialValueFineGrid = 0.
44	!
45	# if ! defined DECAL_FEEDBACK
46	CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
47	! near boundary update:
48	! CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
49	# else
50	CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
51	! near boundary update:
52	! CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
53	# endif
54	!
55	Agrif_UseSpecialValueInUpdate = .FALSE.
56	!
57	#endif
58	!
59	END SUBROUTINE Agrif_Update_Tra
60
61	SUBROUTINE Agrif_Update_Dyn( )
62	!!---------------------------------------------
63	!! * ROUTINE Agrif_Update_Dyn *
64	!!---------------------------------------------
65	!
66	IF (Agrif_Root()) RETURN
67	!
68	#if defined TWO_WAY
69	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update momentum from grid Number',Agrif_Fixed()
70
71	Agrif_UseSpecialValueInUpdate = .FALSE.
72	Agrif_SpecialValueFineGrid = 0.
73	!
74	# if ! defined DECAL_FEEDBACK
75	CALL Agrif_Update_Variable(un_update_id,procname = updateU)
76	CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
77	! near boundary update:
78	! CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
79	! CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
80	# else
81	CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
82	CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
83	! near boundary update:
84	! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
85	! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
86	# endif
87
88	# if ! defined DECAL_FEEDBACK
89	CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
90	CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)
91	# else
92	CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
93	CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)
94	# endif
95	!
96	# if ! defined DECAL_FEEDBACK
97	! Account for updated thicknesses at boundary edges
98	IF (.NOT.ln_linssh) THEN
99	CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
100	CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
101	ENDIF
102	# endif
103	!
104	IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
105	! Update time integrated transports
106	# if ! defined DECAL_FEEDBACK
107	CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
108	CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
109	# else
110	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
111	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
112	# endif
113	END IF
114	#endif
115	!
116	END SUBROUTINE Agrif_Update_Dyn
117
118	SUBROUTINE Agrif_Update_ssh( )
119	!!---------------------------------------------
120	!! * ROUTINE Agrif_Update_ssh *
121	!!---------------------------------------------
122	!
123	IF (Agrif_Root()) RETURN
124	!
125	#if defined TWO_WAY
126	!
127	Agrif_UseSpecialValueInUpdate = .TRUE.
128	Agrif_SpecialValueFineGrid = 0.
129	# if ! defined DECAL_FEEDBACK
130	CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
131	# else
132	CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
133	# endif
134	!
135	Agrif_UseSpecialValueInUpdate = .FALSE.
136	!
137	# if defined DECAL_FEEDBACK && defined VOL_REFLUX
138	IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
139	! Refluxing on ssh:
140	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
141	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
142	END IF
143	# endif
144	!
145	#endif
146	!
147	END SUBROUTINE Agrif_Update_ssh
148
149	# if defined key_zdftke
150
151	SUBROUTINE Agrif_Update_Tke( kt )
152	!!---------------------------------------------
153	!! * ROUTINE Agrif_Update_Tke *
154	!!---------------------------------------------
155	!!
156	INTEGER, INTENT(in) :: kt
157	!
158	IF (Agrif_Root()) RETURN
159	!
160	IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
161	# if defined TWO_WAY
162
163	Agrif_UseSpecialValueInUpdate = .TRUE.
164	Agrif_SpecialValueFineGrid = 0.
165
166	CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN )
167	CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
168	CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
169
170	Agrif_UseSpecialValueInUpdate = .FALSE.
171
172	# endif
173
174	END SUBROUTINE Agrif_Update_Tke
175
176	# endif /* key_zdftke */
177
178	SUBROUTINE Agrif_Update_vvl( )
179	!!---------------------------------------------
180	!! * ROUTINE Agrif_Update_vvl *
181	!!---------------------------------------------
182	!
183	IF (Agrif_Root()) RETURN
184	!
185	#if defined TWO_WAY
186	!
187	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
188	!
189	Agrif_UseSpecialValueInUpdate = .TRUE.
190	Agrif_SpecialValueFineGrid = 0.
191	!
192	! No interface separation here, update vertical grid at T points
193	! everywhere over the overlapping regions (one account for refluxing in that case):
194	CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)
195	!
196	Agrif_UseSpecialValueInUpdate = .FALSE.
197	!
198	CALL Agrif_ChildGrid_To_ParentGrid()
199	CALL dom_vvl_update_UVF
200	CALL Agrif_ParentGrid_To_ChildGrid()
201	!
202	#endif
203	!
204	END SUBROUTINE Agrif_Update_vvl
205
206	SUBROUTINE dom_vvl_update_UVF
207	!!---------------------------------------------
208	!! * ROUTINE dom_vvl_update_UVF *
209	!!---------------------------------------------
210	!!
211	INTEGER :: jk
212	REAL(wp):: zcoef
213	!!---------------------------------------------
214
215	IF (lwp.AND.lk_agrif_debug) Write(,) 'Finalize e3 on grid Number', &
216	& Agrif_Fixed(), 'Step', Agrif_Nb_Step()
217
218	! Save "old" scale factor (prior update) for subsequent asselin correction
219	! of prognostic variables
220	! -----------------------
221	!
222	e3u_a(:,:,:) = e3u_n(:,:,:)
223	e3v_a(:,:,:) = e3v_n(:,:,:)
224	! ua(:,:,:) = e3u_b(:,:,:)
225	! va(:,:,:) = e3v_b(:,:,:)
226	hu_a(:,:) = hu_n(:,:)
227	hv_a(:,:) = hv_n(:,:)
228
229	! 1) NOW fields
230	!--------------
231
232	! Vertical scale factor interpolations
233	! ------------------------------------
234	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) , 'U' )
235	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) , 'V' )
236	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) , 'F' )
237
238	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
239	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
240
241	! Update total depths:
242	! --------------------
243	hu_n(:,:) = 0._wp ! Ocean depth at U-points
244	hv_n(:,:) = 0._wp ! Ocean depth at V-points
245	DO jk = 1, jpkm1
246	hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
247	hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
248	END DO
249	! ! Inverse of the local depth
250	r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
251	r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
252
253
254	! 2) BEFORE fields:
255	!------------------
256	! IF ( (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
257	! & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts &
258	! & .AND.(.NOT.ln_bt_fw)))) THEN
259	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
260	!
261	! Vertical scale factor interpolations
262	! ------------------------------------
263	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
264	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
265
266	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
267	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
268
269	! Update total depths:
270	! --------------------
271	hu_b(:,:) = 0._wp ! Ocean depth at U-points
272	hv_b(:,:) = 0._wp ! Ocean depth at V-points
273	DO jk = 1, jpkm1
274	hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
275	hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
276	END DO
277	! ! Inverse of the local depth
278	r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
279	r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
280	ENDIF
281	!
282	END SUBROUTINE dom_vvl_update_UVF
283
284	SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
285	!!---------------------------------------------
286	!! * ROUTINE updateT *
287	!!---------------------------------------------
288	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
289	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
290	LOGICAL, INTENT(in) :: before
291	!!
292	INTEGER :: ji,jj,jk,jn
293	REAL(wp) :: ztb, ztnu, ztno
294	!!---------------------------------------------
295	!
296	IF (before) THEN
297	DO jn = n1,n2
298	DO jk=k1,k2
299	DO jj=j1,j2
300	DO ji=i1,i2
301	!> jc tmp
302	tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
303	! tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk)
304	!< jc tmp
305	END DO
306	END DO
307	END DO
308	END DO
309	ELSE
310	!> jc tmp
311	DO jn = n1,n2
312	tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
313	& * tmask(i1:i2,j1:j2,k1:k2)
314	ENDDO
315	!< jc tmp
316	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
317	! Add asselin part
318	DO jn = n1,n2
319	DO jk=k1,k2
320	DO jj=j1,j2
321	DO ji=i1,i2
322	IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
323	ztb = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
324	ztnu = tabres(ji,jj,jk,jn)
325	ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
326	tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) &
327	& * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
328	ENDIF
329	ENDDO
330	ENDDO
331	ENDDO
332	ENDDO
333	ENDIF
334	DO jn = n1,n2
335	DO jk=k1,k2
336	DO jj=j1,j2
337	DO ji=i1,i2
338	IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
339	tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
340	END IF
341	END DO
342	END DO
343	END DO
344	END DO
345	!
346	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
347	tsb(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
348	ENDIF
349	!
350	ENDIF
351	!
352	END SUBROUTINE updateTS
353
354
355	SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )
356	!!---------------------------------------------
357	!! * ROUTINE updateu *
358	!!---------------------------------------------
359	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
360	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
361	LOGICAL , INTENT(in ) :: before
362	!
363	INTEGER :: ji, jj, jk
364	REAL(wp) :: zrhoy, zub, zunu, zuno
365	!!---------------------------------------------
366	!
367	IF( before ) THEN
368	zrhoy = Agrif_Rhoy()
369	DO jk = k1, k2
370	tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
371	END DO
372	ELSE
373	DO jk=k1,k2
374	DO jj=j1,j2
375	DO ji=i1,i2
376	tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj)
377	!
378	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
379	zub = ub(ji,jj,jk) * e3u_b(ji,jj,jk) ! fse3t_b prior update should be used
380	zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
381	zunu = tabres(ji,jj,jk)
382	ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &
383	& * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
384	ENDIF
385	!
386	un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
387	END DO
388	END DO
389	END DO
390	!
391	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
392	ub(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
393	ENDIF
394	!
395	ENDIF
396	!
397	END SUBROUTINE updateu
398
399	SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
400	!!---------------------------------------------
401	!! * ROUTINE correct_u_bdy *
402	!!---------------------------------------------
403	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
404	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
405	LOGICAL , INTENT(in ) :: before
406	INTEGER , INTENT(in) :: nb, ndir
407	!!
408	LOGICAL :: western_side, eastern_side
409	!
410	INTEGER :: jj, jk
411	REAL(wp) :: zcor
412	!!---------------------------------------------
413	!
414	IF( .NOT.before ) THEN
415	!
416	western_side = (nb == 1).AND.(ndir == 1)
417	eastern_side = (nb == 1).AND.(ndir == 2)
418	!
419	IF (western_side) THEN
420	DO jj=j1,j2
421	zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
422	un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
423	DO jk=1,jpkm1
424	un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
425	END DO
426	END DO
427	ENDIF
428	!
429	IF (eastern_side) THEN
430	DO jj=j1,j2
431	zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
432	un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
433	DO jk=1,jpkm1
434	un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
435	END DO
436	END DO
437	ENDIF
438	!
439	ENDIF
440	!
441	END SUBROUTINE correct_u_bdy
442
443
444	SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before)
445	!!---------------------------------------------
446	!! * ROUTINE updatev *
447	!!---------------------------------------------
448	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
449	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
450	LOGICAL , INTENT(in ) :: before
451	!
452	INTEGER :: ji, jj, jk
453	REAL(wp) :: zrhox, zvb, zvnu, zvno
454	!!---------------------------------------------
455	!
456	IF (before) THEN
457	zrhox = Agrif_Rhox()
458	DO jk=k1,k2
459	DO jj=j1,j2
460	DO ji=i1,i2
461	tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
462	END DO
463	END DO
464	END DO
465	ELSE
466	DO jk=k1,k2
467	DO jj=j1,j2
468	DO ji=i1,i2
469	tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj)
470	!
471	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
472	zvb = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
473	zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
474	zvnu = tabres(ji,jj,jk)
475	vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &
476	& * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
477	ENDIF
478	!
479	vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
480	END DO
481	END DO
482	END DO
483	!
484	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
485	vb(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
486	ENDIF
487	!
488	ENDIF
489	!
490	END SUBROUTINE updatev
491
492	SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
493	!!---------------------------------------------
494	!! * ROUTINE correct_u_bdy *
495	!!---------------------------------------------
496	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
497	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
498	LOGICAL , INTENT(in ) :: before
499	INTEGER , INTENT(in) :: nb, ndir
500	!!
501	LOGICAL :: southern_side, northern_side
502	!
503	INTEGER :: ji, jk
504	REAL(wp) :: zcor
505	!!---------------------------------------------
506	!
507	IF( .NOT.before ) THEN
508	!
509	southern_side = (nb == 2).AND.(ndir == 1)
510	northern_side = (nb == 2).AND.(ndir == 2)
511	!
512	IF (southern_side) THEN
513	DO ji=i1,i2
514	zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
515	vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
516	DO jk=1,jpkm1
517	vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
518	END DO
519	END DO
520	ENDIF
521	!
522	IF (northern_side) THEN
523	DO ji=i1,i2
524	zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
525	vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
526	DO jk=1,jpkm1
527	vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
528	END DO
529	END DO
530	ENDIF
531	!
532	ENDIF
533	!
534	END SUBROUTINE correct_v_bdy
535
536
537	SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
538	!!---------------------------------------------
539	!! * ROUTINE updateu2d *
540	!!---------------------------------------------
541	INTEGER, INTENT(in) :: i1, i2, j1, j2
542	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
543	LOGICAL, INTENT(in) :: before
544	!!
545	INTEGER :: ji, jj, jk
546	REAL(wp) :: zrhoy
547	REAL(wp) :: zcorr
548	!!---------------------------------------------
549	!
550	IF (before) THEN
551	zrhoy = Agrif_Rhoy()
552	DO jj=j1,j2
553	DO ji=i1,i2
554	tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
555	END DO
556	END DO
557	ELSE
558	DO jj=j1,j2
559	DO ji=i1,i2
560	tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj)
561	!
562	! Update "now" 3d velocities:
563	spgu(ji,jj) = 0._wp
564	DO jk=1,jpkm1
565	spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
566	END DO
567	!
568	zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
569	DO jk=1,jpkm1
570	un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)
571	END DO
572	!
573	! Update barotropic velocities:
574	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
575	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
576	zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
577	ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
578	END IF
579	ENDIF
580	un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
581	!
582	! Correct "before" velocities to hold correct bt component:
583	spgu(ji,jj) = 0.e0
584	DO jk=1,jpkm1
585	spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
586	END DO
587	!
588	zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
589	DO jk=1,jpkm1
590	ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)
591	END DO
592	!
593	END DO
594	END DO
595	!
596	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
597	ub_b(i1:i2,j1:j2) = un_b(i1:i2,j1:j2)
598	ENDIF
599	ENDIF
600	!
601	END SUBROUTINE updateu2d
602
603
604	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
605	!!---------------------------------------------
606	!! * ROUTINE updatev2d *
607	!!---------------------------------------------
608	INTEGER, INTENT(in) :: i1, i2, j1, j2
609	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
610	LOGICAL, INTENT(in) :: before
611	!!
612	INTEGER :: ji, jj, jk
613	REAL(wp) :: zrhox
614	REAL(wp) :: zcorr
615	!!---------------------------------------------
616	!
617	IF (before) THEN
618	zrhox = Agrif_Rhox()
619	DO jj=j1,j2
620	DO ji=i1,i2
621	tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj)
622	END DO
623	END DO
624	ELSE
625	DO jj=j1,j2
626	DO ji=i1,i2
627	tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj)
628	!
629	! Update "now" 3d velocities:
630	spgv(ji,jj) = 0.e0
631	DO jk=1,jpkm1
632	spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
633	END DO
634	!
635	zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
636	DO jk=1,jpkm1
637	vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)
638	END DO
639	!
640	! Update barotropic velocities:
641	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
642	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
643	zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
644	vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
645	END IF
646	ENDIF
647	vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
648	!
649	! Correct "before" velocities to hold correct bt component:
650	spgv(ji,jj) = 0.e0
651	DO jk=1,jpkm1
652	spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
653	END DO
654	!
655	zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
656	DO jk=1,jpkm1
657	vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)
658	END DO
659	!
660	END DO
661	END DO
662	!
663	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
664	vb_b(i1:i2,j1:j2) = vn_b(i1:i2,j1:j2)
665	ENDIF
666	!
667	ENDIF
668	!
669	END SUBROUTINE updatev2d
670
671
672	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
673	!!---------------------------------------------
674	!! * ROUTINE updateSSH *
675	!!---------------------------------------------
676	INTEGER, INTENT(in) :: i1, i2, j1, j2
677	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
678	LOGICAL, INTENT(in) :: before
679	!!
680	INTEGER :: ji, jj
681	!!---------------------------------------------
682	!
683	IF (before) THEN
684	DO jj=j1,j2
685	DO ji=i1,i2
686	tabres(ji,jj) = sshn(ji,jj)
687	END DO
688	END DO
689	ELSE
690	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
691	DO jj=j1,j2
692	DO ji=i1,i2
693	sshb(ji,jj) = sshb(ji,jj) &
694	& + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
695	END DO
696	END DO
697	ENDIF
698	!
699	DO jj=j1,j2
700	DO ji=i1,i2
701	sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
702	END DO
703	END DO
704	!
705	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
706	sshb(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
707	ENDIF
708	!
709
710	ENDIF
711	!
712	END SUBROUTINE updateSSH
713
714
715	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
716	!!---------------------------------------------
717	!! * ROUTINE updateub2b *
718	!!---------------------------------------------
719	INTEGER, INTENT(in) :: i1, i2, j1, j2
720	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
721	LOGICAL, INTENT(in) :: before
722	!!
723	INTEGER :: ji, jj
724	REAL(wp) :: zrhoy, za1, zcor
725	!!---------------------------------------------
726	!
727	IF (before) THEN
728	zrhoy = Agrif_Rhoy()
729	DO jj=j1,j2
730	DO ji=i1,i2
731	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
732	END DO
733	END DO
734	tabres = zrhoy * tabres
735	ELSE
736	!
737	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
738	!
739	za1 = 1._wp / REAL(Agrif_rhot(), wp)
740	DO jj=j1,j2
741	DO ji=i1,i2
742	zcor=tabres(ji,jj) - ub2_b(ji,jj)
743	! Update time integrated fluxes also in case of multiply nested grids:
744	ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor
745	! Update corrective fluxes:
746	un_bf(ji,jj) = un_bf(ji,jj) + zcor
747	! Update half step back fluxes:
748	ub2_b(ji,jj) = tabres(ji,jj)
749	END DO
750	END DO
751	ENDIF
752	!
753	END SUBROUTINE updateub2b
754
755	SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
756	!!---------------------------------------------
757	!! * ROUTINE reflux_sshu *
758	!!---------------------------------------------
759	INTEGER, INTENT(in) :: i1, i2, j1, j2
760	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
761	LOGICAL, INTENT(in) :: before
762	INTEGER, INTENT(in) :: nb, ndir
763	!!
764	LOGICAL :: western_side, eastern_side
765	INTEGER :: ji, jj
766	REAL(wp) :: zrhoy, za1, zcor
767	!!---------------------------------------------
768	!
769	IF (before) THEN
770	zrhoy = Agrif_Rhoy()
771	DO jj=j1,j2
772	DO ji=i1,i2
773	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
774	END DO
775	END DO
776	tabres = zrhoy * tabres
777	ELSE
778	!
779	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
780	!
781	western_side = (nb == 1).AND.(ndir == 1)
782	eastern_side = (nb == 1).AND.(ndir == 2)
783	!
784	IF (western_side) THEN
785	DO jj=j1,j2
786	zcor = rdt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
787	sshn(i1 ,jj) = sshn(i1 ,jj) + zcor
788	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1 ,jj) = sshb(i1 ,jj) + atfp * zcor
789	END DO
790	ENDIF
791	IF (eastern_side) THEN
792	DO jj=j1,j2
793	zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
794	sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
795	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
796	END DO
797	ENDIF
798	!
799	ENDIF
800	!
801	END SUBROUTINE reflux_sshu
802
803	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
804	!!---------------------------------------------
805	!! * ROUTINE updatevb2b *
806	!!---------------------------------------------
807	INTEGER, INTENT(in) :: i1, i2, j1, j2
808	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
809	LOGICAL, INTENT(in) :: before
810	!!
811	INTEGER :: ji, jj
812	REAL(wp) :: zrhox, za1, zcor
813	!!---------------------------------------------
814	!
815	IF (before) THEN
816	zrhox = Agrif_Rhox()
817	DO jj=j1,j2
818	DO ji=i1,i2
819	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
820	END DO
821	END DO
822	tabres = zrhox * tabres
823	ELSE
824	!
825	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
826	!
827	za1 = 1._wp / REAL(Agrif_rhot(), wp)
828	DO jj=j1,j2
829	DO ji=i1,i2
830	zcor=tabres(ji,jj) - vb2_b(ji,jj)
831	! Update time integrated fluxes also in case of multiply nested grids:
832	vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor
833	! Update corrective fluxes:
834	vn_bf(ji,jj) = vn_bf(ji,jj) + zcor
835	! Update half step back fluxes:
836	vb2_b(ji,jj) = tabres(ji,jj)
837	END DO
838	END DO
839	ENDIF
840	!
841	END SUBROUTINE updatevb2b
842
843	SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
844	!!---------------------------------------------
845	!! * ROUTINE reflux_sshv *
846	!!---------------------------------------------
847	INTEGER, INTENT(in) :: i1, i2, j1, j2
848	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
849	LOGICAL, INTENT(in) :: before
850	INTEGER, INTENT(in) :: nb, ndir
851	!!
852	LOGICAL :: southern_side, northern_side
853	INTEGER :: ji, jj
854	REAL(wp) :: zrhox, za1, zcor
855	!!---------------------------------------------
856	!
857	IF (before) THEN
858	zrhox = Agrif_Rhox()
859	DO jj=j1,j2
860	DO ji=i1,i2
861	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
862	END DO
863	END DO
864	tabres = zrhox * tabres
865	ELSE
866	!
867	tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
868	!
869	southern_side = (nb == 2).AND.(ndir == 1)
870	northern_side = (nb == 2).AND.(ndir == 2)
871	!
872	IF (southern_side) THEN
873	DO ji=i1,i2
874	zcor = rdt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1))
875	sshn(ji,j1 ) = sshn(ji,j1 ) + zcor
876	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1 ) = sshb(ji,j1) + atfp * zcor
877	END DO
878	ENDIF
879	IF (northern_side) THEN
880	DO ji=i1,i2
881	zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2))
882	sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
883	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
884	END DO
885	ENDIF
886	!
887	ENDIF
888	!
889	END SUBROUTINE reflux_sshv
890
891	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
892	! currently not used
893	!!---------------------------------------------
894	!! * ROUTINE updateT *
895	!!---------------------------------------------
896	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
897	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
898	LOGICAL, iNTENT(in) :: before
899	!
900	INTEGER :: ji,jj,jk
901	REAL(wp) :: ztemp
902	!!---------------------------------------------
903
904	IF (before) THEN
905	DO jk=k1,k2
906	DO jj=j1,j2
907	DO ji=i1,i2
908	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
909	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
910	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
911	END DO
912	END DO
913	END DO
914	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
915	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
916	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
917	ELSE
918	DO jk=k1,k2
919	DO jj=j1,j2
920	DO ji=i1,i2
921	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
922	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
923	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
924	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
925	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
926	print *,'CORR = ',ztemp-1.
927	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
928	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
929	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
930	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
931	END IF
932	END DO
933	END DO
934	END DO
935	ENDIF
936	!
937	END SUBROUTINE update_scales
938
939	# if defined key_zdftke
940
941	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
942	!!---------------------------------------------
943	!! * ROUTINE updateen *
944	!!---------------------------------------------
945	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
946	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
947	LOGICAL, INTENT(in) :: before
948	!!---------------------------------------------
949	!
950	IF (before) THEN
951	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
952	ELSE
953	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
954	ENDIF
955	!
956	END SUBROUTINE updateEN
957
958
959	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
960	!!---------------------------------------------
961	!! * ROUTINE updateavt *
962	!!---------------------------------------------
963	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
964	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
965	LOGICAL, INTENT(in) :: before
966	!!---------------------------------------------
967	!
968	IF (before) THEN
969	ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
970	ELSE
971	avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
972	ENDIF
973	!
974	END SUBROUTINE updateAVT
975
976
977	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
978	!!---------------------------------------------
979	!! * ROUTINE updateavm *
980	!!---------------------------------------------
981	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
982	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
983	LOGICAL, INTENT(in) :: before
984	!!---------------------------------------------
985	!
986	IF (before) THEN
987	ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
988	ELSE
989	avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
990	ENDIF
991	!
992	END SUBROUTINE updateAVM
993
994	# endif /* key_zdftke */
995
996	SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
997	!!---------------------------------------------
998	!! * ROUTINE updatee3t *
999	!!---------------------------------------------
1000	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1001	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1002	LOGICAL, INTENT(in) :: before
1003	!
1004	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1005	INTEGER :: ji,jj,jk
1006	REAL(wp) :: zcoef
1007	!!---------------------------------------------
1008	!
1009	IF (.NOT.before) THEN
1010	!
1011	ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))
1012	!
1013	! Update e3t from ssh (z* case only)
1014	DO jk = 1, jpkm1
1015	DO jj=j1,j2
1016	DO ji=i1,i2
1017	ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
1018	& *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1019	END DO
1020	END DO
1021	END DO
1022	!
1023	! 1) Updates at BEFORE time step:
1024	! -------------------------------
1025	!
1026	! Save "old" scale factor (prior update) for subsequent asselin correction
1027	! of prognostic variables
1028	e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1029
1030	! One should also save e3t_b, but lacking of workspace...
1031	! hdivn(i1:i2,j1:j2,1:jpkm1) = e3t_b(i1:i2,j1:j2,1:jpkm1)
1032
1033	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1034	DO jk = 1, jpkm1
1035	DO jj=j1,j2
1036	DO ji=i1,i2
1037	e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) &
1038	& + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1039	END DO
1040	END DO
1041	END DO
1042	!
1043	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)
1044	gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1045	gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1046	!
1047	DO jk = 2, jpk
1048	DO jj = j1,j2
1049	DO ji = i1,i2
1050	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1051	e3w_b(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &
1052	& ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) &
1053	& + 0.5_wp * tmask(ji,jj,jk) * &
1054	& ( e3t_b(ji,jj,jk ) - e3t_0(ji,jj,jk ) )
1055	gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1056	gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &
1057	& + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk))
1058	END DO
1059	END DO
1060	END DO
1061	!
1062	ENDIF
1063	!
1064	! 2) Updates at NOW time step:
1065	! ----------------------------
1066	!
1067	! Update vertical scale factor at T-points:
1068	e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1069	!
1070	! Update total depth:
1071	ht_n(i1:i2,j1:j2) = 0._wp
1072	DO jk = 1, jpkm1
1073	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)
1074	END DO
1075	!
1076	! Update vertical scale factor at W-points and depths:
1077	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)
1078	gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1079	gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1080	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
1081	!
1082	DO jk = 2, jpk
1083	DO jj = j1,j2
1084	DO ji = i1,i2
1085	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1086	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) ) &
1087	& + 0.5_wp * tmask(ji,jj,jk) * ( e3t_n(ji,jj,jk ) - e3t_0(ji,jj,jk ) )
1088	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1089	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &
1090	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
1091	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
1092	END DO
1093	END DO
1094	END DO
1095	!
1096	IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1097	e3t_b (i1:i2,j1:j2,1:jpk) = e3t_n (i1:i2,j1:j2,1:jpk)
1098	e3w_b (i1:i2,j1:j2,1:jpk) = e3w_n (i1:i2,j1:j2,1:jpk)
1099	gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1100	gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1101	ENDIF
1102	!
1103	DEALLOCATE(ptab)
1104	ENDIF
1105	!
1106	END SUBROUTINE updatee3t
1107
1108	#else
1109	CONTAINS
1110	SUBROUTINE agrif_opa_update_empty
1111	!!---------------------------------------------
1112	!! * ROUTINE agrif_opa_update_empty *
1113	!!---------------------------------------------
1114	WRITE(,) 'agrif_opa_update : You should not have seen this print! error?'
1115	END SUBROUTINE agrif_opa_update_empty
1116	#endif
1117	END MODULE agrif_opa_update
1118

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

Download in other formats: