New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_oce_update.F90 in NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/NST – NEMO

source: NEMO/branches/2020/dev_r12973_AGRIF_CMEMS/src/NST/agrif_oce_update.F90 @ 13026

Last change on this file since 13026 was 13026, checked in by rblod, 3 years ago

AGRIF with northfold and perio, see ticket #2129

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