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

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

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

Download in other formats: