source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90 @ 11827

Last change on this file since 11827 was 11827, checked in by jchanut, 11 months ago

#2222, corrections to interpolation within sponge

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