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

Last change on this file since 11760 was 11760, checked in by jchanut, 17 months ago

#2222, mistake in previous commit

  • 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! jc_alt      Agrif_UseSpecialValueInUpdate = .FALSE.
52      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         AGRIF_SpecialValue = -999._wp
298         DO jn = n1,n2-1
299            DO jk=k1,k2
300               DO jj=j1,j2
301                  DO ji=i1,i2
302                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
303                                         &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp
304!jc_alt                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk)
305                  END DO
306               END DO
307            END DO
308         END DO
309         DO jk=k1,k2
310            DO jj=j1,j2
311               DO ji=i1,i2
312                  tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
313                                      &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp
314! jc_alt                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
315               END DO
316            END DO
317         END DO
318      ELSE
319         tabres_child(:,:,:,:) = 0._wp
320         AGRIF_SpecialValue = 0._wp
321         DO jj=j1,j2
322            DO ji=i1,i2
323               N_in = 0
324               DO jk=k1,k2 !k2 = jpk of child grid
325                  IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT
326! jc_alt                  IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT
327                  N_in = N_in + 1
328                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
329                  h_in(N_in) = tabres(ji,jj,jk,n2)
330               ENDDO
331               N_out = 0
332               DO jk=1,jpk ! jpk of parent grid
333                  IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
334                  N_out = N_out + 1
335                  h_out(N_out) = e3t_n(ji,jj,jk) 
336               ENDDO
337               IF (N_in*N_out > 0) THEN !Remove this?
338                  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)
339               ENDIF
340            ENDDO
341         ENDDO
342
343         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
344            ! Add asselin part
345            DO jn = 1,jpts
346               DO jk = 1, jpkm1
347                  DO jj = j1, j2
348                     DO ji = i1, i2
349                        IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
350                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
351                           ztnu = tabres_child(ji,jj,jk,jn) * e3t_n(ji,jj,jk)
352                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
353                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
354                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
355                        ENDIF
356                     END DO
357                  END DO
358               END DO
359            END DO
360         ENDIF
361         DO jn = 1,jpts
362            DO jk = 1, jpkm1
363               DO jj = j1, j2
364                  DO ji = i1, i2
365                     IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
366                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn)
367                     END IF
368                  END DO
369               END DO
370            END DO
371         END DO
372         !
373         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
374            tsb(i1:i2,j1:j2,1:jpkm1,1:jpts)  = tsn(i1:i2,j1:j2,1:jpkm1,1:jpts)
375         ENDIF
376      ENDIF
377      !
378   END SUBROUTINE updateTS
379
380# else
381
382   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
383      !!---------------------------------------------
384      !!           *** ROUTINE updateT ***
385      !!---------------------------------------------
386      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
387      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
388      LOGICAL, INTENT(in) :: before
389      !!
390      INTEGER :: ji,jj,jk,jn
391      REAL(wp) :: ztb, ztnu, ztno
392      !!---------------------------------------------
393      !
394      IF (before) THEN
395         DO jn = 1,jpts
396            DO jk=k1,k2
397               DO jj=j1,j2
398                  DO ji=i1,i2
399!> jc tmp
400                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
401!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
402!< jc tmp
403                  END DO
404               END DO
405            END DO
406         END DO
407      ELSE
408!> jc tmp
409         DO jn = 1,jpts
410            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
411                                         & * tmask(i1:i2,j1:j2,k1:k2)
412         ENDDO
413!< jc tmp
414         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
415            ! Add asselin part
416            DO jn = 1,jpts
417               DO jk = k1, k2
418                  DO jj = j1, j2
419                     DO ji = i1, i2
420                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
421                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
422                           ztnu = tabres(ji,jj,jk,jn)
423                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
424                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
425                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
426                        ENDIF
427                     END DO
428                  END DO
429               END DO
430            END DO
431         ENDIF
432         DO jn = 1,jpts
433            DO jk=k1,k2
434               DO jj=j1,j2
435                  DO ji=i1,i2
436                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
437                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
438                     END IF
439                  END DO
440               END DO
441            END DO
442         END DO
443         !
444         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
445            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
446         ENDIF
447         !
448      ENDIF
449      !
450   END SUBROUTINE updateTS
451
452# endif
453
454# if defined key_vertical
455
456   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
457      !!---------------------------------------------
458      !!           *** ROUTINE updateu ***
459      !!---------------------------------------------
460      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
461      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
462      LOGICAL                                     , INTENT(in   ) :: before
463      !
464      INTEGER ::   ji, jj, jk
465      REAL(wp)::   zrhoy, zub, zunu, zuno
466! VERTICAL REFINEMENT BEGIN
467      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
468      REAL(wp) :: h_in(k1:k2)
469      REAL(wp) :: h_out(1:jpk)
470      INTEGER  :: N_in, N_out
471      REAL(wp) :: h_diff, excess, thick
472      REAL(wp) :: tabin(k1:k2)
473! VERTICAL REFINEMENT END
474      !!---------------------------------------------
475      !
476      IF( before ) THEN
477         zrhoy = Agrif_Rhoy()
478         AGRIF_SpecialValue = -999._wp
479         DO jk=k1,k2
480            DO jj=j1,j2
481               DO ji=i1,i2
482                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
483                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp
484! jc_alt                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) 
485                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
486                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp
487! jc_alt                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
488               END DO
489            END DO
490         END DO
491      ELSE
492         tabres_child(:,:,:) = 0.
493         AGRIF_SpecialValue = 0._wp
494         DO jj=j1,j2
495            DO ji=i1,i2
496               N_in = 0
497               h_in(:) = 0._wp
498               tabin(:) = 0._wp
499               DO jk=k1,k2 !k2=jpk of child grid
500                  IF( tabres(ji,jj,jk,2) < -900._wp) EXIT
501! jc_alt                  IF( tabres(ji,jj,jk,2) == 0.) EXIT
502                  N_in = N_in + 1
503                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
504                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
505               ENDDO
506               N_out = 0
507               DO jk=1,jpk
508                  IF (umask(ji,jj,jk) == 0) EXIT
509                  N_out = N_out + 1
510                  h_out(N_out) = e3u_n(ji,jj,jk)
511               ENDDO
512               IF (N_in * N_out > 0) THEN
513                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
514                  IF (h_diff < -1.e-4) THEN
515!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
516!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
517                     excess = 0._wp
518                     DO jk=N_in,1,-1
519                        thick = MIN(-1*h_diff, h_in(jk))
520                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
521                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
522                        h_diff = h_diff + thick
523                        IF ( h_diff == 0) THEN
524                           N_in = jk
525                           h_in(jk) = h_in(jk) - thick
526                           EXIT
527                        ENDIF
528                     ENDDO
529                  ENDIF
530                  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)
531                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
532               ENDIF
533            ENDDO
534         ENDDO
535         !
536         DO jk=1,jpk
537            DO jj=j1,j2
538               DO ji=i1,i2
539                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
540                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
541                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
542                     zunu = tabres_child(ji,jj,jk) * e3u_n(ji,jj,jk)
543                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &     
544                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
545                  ENDIF
546                  !
547                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
548               END DO
549            END DO
550         END DO
551         !
552         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
553            ub(i1:i2,j1:j2,1:jpkm1)  = un(i1:i2,j1:j2,1:jpkm1)
554         ENDIF
555         !
556      ENDIF
557      !
558   END SUBROUTINE updateu
559
560#else
561
562   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
563      !!---------------------------------------------
564      !!           *** ROUTINE updateu ***
565      !!---------------------------------------------
566      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
567      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
568      LOGICAL                                     , INTENT(in   ) :: before
569      !
570      INTEGER  :: ji, jj, jk
571      REAL(wp) :: zrhoy, zub, zunu, zuno
572      !!---------------------------------------------
573      !
574      IF( before ) THEN
575         zrhoy = Agrif_Rhoy()
576         DO jk = k1, k2
577            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)
578         END DO
579      ELSE
580         DO jk=k1,k2
581            DO jj=j1,j2
582               DO ji=i1,i2
583                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 
584                  !
585                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
586                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
587                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
588                     zunu = tabres(ji,jj,jk,1)
589                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &     
590                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
591                  ENDIF
592                  !
593                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
594               END DO
595            END DO
596         END DO
597         !
598         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
599            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
600         ENDIF
601         !
602      ENDIF
603      !
604   END SUBROUTINE updateu
605
606# endif
607
608   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
609      !!---------------------------------------------
610      !!           *** ROUTINE correct_u_bdy ***
611      !!---------------------------------------------
612      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
613      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
614      LOGICAL                                     , INTENT(in   ) :: before
615      INTEGER                                     , INTENT(in)    :: nb, ndir
616      !!
617      LOGICAL :: western_side, eastern_side 
618      !
619      INTEGER  :: jj, jk
620      REAL(wp) :: zcor
621      !!---------------------------------------------
622      !
623      IF( .NOT.before ) THEN
624         !
625         western_side  = (nb == 1).AND.(ndir == 1)
626         eastern_side  = (nb == 1).AND.(ndir == 2)
627         !
628         IF (western_side) THEN
629            DO jj=j1,j2
630               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
631               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
632               DO jk=1,jpkm1
633                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
634               END DO
635            END DO
636         ENDIF
637         !
638         IF (eastern_side) THEN
639            DO jj=j1,j2
640               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
641               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
642               DO jk=1,jpkm1
643                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
644               END DO
645            END DO
646         ENDIF
647         !
648      ENDIF
649      !
650   END SUBROUTINE correct_u_bdy
651
652# if  defined key_vertical
653
654   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
655      !!---------------------------------------------
656      !!           *** ROUTINE updatev ***
657      !!---------------------------------------------
658      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
659      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
660      LOGICAL                                     , INTENT(in   ) :: before
661      !
662      INTEGER  ::   ji, jj, jk
663      REAL(wp) ::   zrhox, zvb, zvnu, zvno
664! VERTICAL REFINEMENT BEGIN
665      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
666      REAL(wp) :: h_in(k1:k2)
667      REAL(wp) :: h_out(1:jpk)
668      INTEGER :: N_in, N_out
669      REAL(wp) :: h_diff, excess, thick
670      REAL(wp) :: tabin(k1:k2)
671! VERTICAL REFINEMENT END
672      !!---------------------------------------------     
673      !
674      IF( before ) THEN
675         zrhox = Agrif_Rhox()
676         AGRIF_SpecialValue = -999._wp
677         DO jk=k1,k2
678            DO jj=j1,j2
679               DO ji=i1,i2
680                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
681                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp
682! jc_alt                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)
683                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) &
684                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp
685! jc_alt                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
686               END DO
687            END DO
688         END DO
689      ELSE
690         tabres_child(:,:,:) = 0.
691         AGRIF_SpecialValue = 0._wp
692         DO jj=j1,j2
693            DO ji=i1,i2
694               N_in = 0
695               DO jk=k1,k2
696                  IF (tabres(ji,jj,jk,2) < -900._wp) EXIT
697! jc_alt                  IF (tabres(ji,jj,jk,2) == 0) EXIT
698                  N_in = N_in + 1
699                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
700                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
701               ENDDO
702               N_out = 0
703               DO jk=1,jpk
704                  IF (vmask(ji,jj,jk) == 0) EXIT
705                  N_out = N_out + 1
706                  h_out(N_out) = e3v_n(ji,jj,jk)
707               ENDDO
708               IF (N_in * N_out > 0) THEN
709                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
710                  IF (h_diff < -1.e-4) then
711!Even if bathy at T points match it's possible for the V points to be deeper in the child grid.
712!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
713                     excess = 0._wp
714                     DO jk=N_in,1,-1
715                        thick = MIN(-1*h_diff, h_in(jk))
716                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
717                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
718                        h_diff = h_diff + thick
719                        IF ( h_diff == 0) THEN
720                           N_in = jk
721                           h_in(jk) = h_in(jk) - thick
722                           EXIT
723                        ENDIF
724                     ENDDO
725                  ENDIF
726                  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)
727                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
728               ENDIF
729            ENDDO
730         ENDDO
731         !
732         DO jk=1,jpkm1
733            DO jj=j1,j2
734               DO ji=i1,i2
735                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
736                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
737                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
738                     zvnu = tabres_child(ji,jj,jk) * e3v_n(ji,jj,jk)
739                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &     
740                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
741                  ENDIF
742                  !
743                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
744               END DO
745            END DO
746         END DO
747         !
748         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
749            vb(i1:i2,j1:j2,1:jpkm1)  = vn(i1:i2,j1:j2,1:jpkm1)
750         ENDIF
751         !
752      ENDIF
753      !
754   END SUBROUTINE updatev
755
756# else
757
758   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
759      !!---------------------------------------------
760      !!           *** ROUTINE updatev ***
761      !!---------------------------------------------
762      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
763      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
764      LOGICAL                                     , INTENT(in   ) :: before
765      !
766      INTEGER  :: ji, jj, jk
767      REAL(wp) :: zrhox, zvb, zvnu, zvno
768      !!---------------------------------------------     
769      !
770      IF (before) THEN
771         zrhox = Agrif_Rhox()
772         DO jk=k1,k2
773            DO jj=j1,j2
774               DO ji=i1,i2
775                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
776               END DO
777            END DO
778         END DO
779      ELSE
780         DO jk=k1,k2
781            DO jj=j1,j2
782               DO ji=i1,i2
783                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
784                  !
785                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
786                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
787                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
788                     zvnu = tabres(ji,jj,jk,1)
789                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &     
790                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
791                  ENDIF
792                  !
793                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
794               END DO
795            END DO
796         END DO
797         !
798         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
799            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
800         ENDIF
801         !
802      ENDIF
803      !
804   END SUBROUTINE updatev
805
806# endif
807
808   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
809      !!---------------------------------------------
810      !!           *** ROUTINE correct_u_bdy ***
811      !!---------------------------------------------
812      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
813      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
814      LOGICAL                                     , INTENT(in   ) :: before
815      INTEGER                                     , INTENT(in)    :: nb, ndir
816      !!
817      LOGICAL :: southern_side, northern_side 
818      !
819      INTEGER  :: ji, jk
820      REAL(wp) :: zcor
821      !!---------------------------------------------
822      !
823      IF( .NOT.before ) THEN
824         !
825         southern_side = (nb == 2).AND.(ndir == 1)
826         northern_side = (nb == 2).AND.(ndir == 2)
827         !
828         IF (southern_side) THEN
829            DO ji=i1,i2
830               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
831               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
832               DO jk=1,jpkm1
833                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
834               END DO
835            END DO
836         ENDIF
837         !
838         IF (northern_side) THEN
839            DO ji=i1,i2
840               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
841               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
842               DO jk=1,jpkm1
843                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
844               END DO
845            END DO
846         ENDIF
847         !
848      ENDIF
849      !
850   END SUBROUTINE correct_v_bdy
851
852
853   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
854      !!----------------------------------------------------------------------
855      !!                      *** ROUTINE updateu2d ***
856      !!----------------------------------------------------------------------
857      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
858      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
859      LOGICAL                         , INTENT(in   ) ::   before
860      !!
861      INTEGER  :: ji, jj, jk
862      REAL(wp) :: zrhoy
863      REAL(wp) :: zcorr
864      !!---------------------------------------------
865      !
866      IF( before ) THEN
867         zrhoy = Agrif_Rhoy()
868         DO jj=j1,j2
869            DO ji=i1,i2
870               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
871            END DO
872         END DO
873      ELSE
874         DO jj=j1,j2
875            DO ji=i1,i2
876               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj) 
877               !   
878               ! Update "now" 3d velocities:
879               spgu(ji,jj) = 0._wp
880               DO jk=1,jpkm1
881                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
882               END DO
883               !
884               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
885               DO jk=1,jpkm1             
886                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
887               END DO
888               !
889               ! Update barotropic velocities:
890               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
891                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
892                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
893                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
894                  END IF
895               ENDIF   
896               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
897               !       
898               ! Correct "before" velocities to hold correct bt component:
899               spgu(ji,jj) = 0.e0
900               DO jk=1,jpkm1
901                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
902               END DO
903               !
904               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
905               DO jk=1,jpkm1             
906                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
907               END DO
908               !
909            END DO
910         END DO
911         !
912         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
913            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
914         ENDIF
915      ENDIF
916      !
917   END SUBROUTINE updateu2d
918
919
920   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
921      !!----------------------------------------------------------------------
922      !!                   *** ROUTINE updatev2d ***
923      !!----------------------------------------------------------------------
924      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
925      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
926      LOGICAL                         , INTENT(in   ) ::   before
927      !
928      INTEGER  :: ji, jj, jk
929      REAL(wp) :: zrhox, zcorr
930      !!----------------------------------------------------------------------
931      !
932      IF( before ) THEN
933         zrhox = Agrif_Rhox()
934         DO jj=j1,j2
935            DO ji=i1,i2
936               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
937            END DO
938         END DO
939      ELSE
940         DO jj=j1,j2
941            DO ji=i1,i2
942               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj) 
943               !   
944               ! Update "now" 3d velocities:
945               spgv(ji,jj) = 0.e0
946               DO jk=1,jpkm1
947                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
948               END DO
949               !
950               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
951               DO jk=1,jpkm1             
952                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
953               END DO
954               !
955               ! Update barotropic velocities:
956               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
957                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
958                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
959                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
960                  END IF
961               ENDIF             
962               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
963               !       
964               ! Correct "before" velocities to hold correct bt component:
965               spgv(ji,jj) = 0.e0
966               DO jk=1,jpkm1
967                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
968               END DO
969               !
970               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
971               DO jk=1,jpkm1             
972                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
973               END DO
974               !
975            END DO
976         END DO
977         !
978         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
979            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
980         ENDIF
981         !
982      ENDIF
983      !
984   END SUBROUTINE updatev2d
985
986
987   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
988      !!----------------------------------------------------------------------
989      !!                   *** ROUTINE updateSSH ***
990      !!----------------------------------------------------------------------
991      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
992      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
993      LOGICAL                         , INTENT(in   ) ::   before
994      !!
995      INTEGER :: ji, jj
996      !!----------------------------------------------------------------------
997      !
998      IF( before ) THEN
999         DO jj=j1,j2
1000            DO ji=i1,i2
1001               tabres(ji,jj) = sshn(ji,jj)
1002            END DO
1003         END DO
1004      ELSE
1005         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
1006            DO jj=j1,j2
1007               DO ji=i1,i2
1008                  sshb(ji,jj) =   sshb(ji,jj) &
1009                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
1010               END DO
1011            END DO
1012         ENDIF
1013         !
1014         DO jj=j1,j2
1015            DO ji=i1,i2
1016               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
1017            END DO
1018         END DO
1019         !
1020         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1021            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
1022         ENDIF
1023         !
1024
1025      ENDIF
1026      !
1027   END SUBROUTINE updateSSH
1028
1029
1030   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1031      !!----------------------------------------------------------------------
1032      !!                      *** ROUTINE updateub2b ***
1033      !!----------------------------------------------------------------------
1034      INTEGER                            , INTENT(in) ::   i1, i2, j1, j2
1035      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
1036      LOGICAL                            , INTENT(in) ::   before
1037      !!
1038      INTEGER :: ji, jj
1039      REAL(wp) :: zrhoy, za1, zcor
1040      !!---------------------------------------------
1041      !
1042      IF (before) THEN
1043         zrhoy = Agrif_Rhoy()
1044         DO jj=j1,j2
1045            DO ji=i1,i2
1046               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1047            END DO
1048         END DO
1049         tabres = zrhoy * tabres
1050      ELSE
1051         !
1052         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1053         !
1054         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1055         DO jj=j1,j2
1056            DO ji=i1,i2
1057               zcor=tabres(ji,jj) - ub2_b(ji,jj)
1058               ! Update time integrated fluxes also in case of multiply nested grids:
1059               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
1060               ! Update corrective fluxes:
1061               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
1062               ! Update half step back fluxes:
1063               ub2_b(ji,jj) = tabres(ji,jj)
1064            END DO
1065         END DO
1066      ENDIF
1067      !
1068   END SUBROUTINE updateub2b
1069
1070   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1071      !!---------------------------------------------
1072      !!          *** ROUTINE reflux_sshu ***
1073      !!---------------------------------------------
1074      INTEGER, INTENT(in) :: i1, i2, j1, j2
1075      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1076      LOGICAL, INTENT(in) :: before
1077      INTEGER, INTENT(in) :: nb, ndir
1078      !!
1079      LOGICAL :: western_side, eastern_side 
1080      INTEGER :: ji, jj
1081      REAL(wp) :: zrhoy, za1, zcor
1082      !!---------------------------------------------
1083      !
1084      IF (before) THEN
1085         zrhoy = Agrif_Rhoy()
1086         DO jj=j1,j2
1087            DO ji=i1,i2
1088               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1089            END DO
1090         END DO
1091         tabres = zrhoy * tabres
1092      ELSE
1093         !
1094         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1095         !
1096         western_side  = (nb == 1).AND.(ndir == 1)
1097         eastern_side  = (nb == 1).AND.(ndir == 2)
1098         !
1099         IF (western_side) THEN
1100            DO jj=j1,j2
1101               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
1102               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
1103               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
1104            END DO
1105         ENDIF
1106         IF (eastern_side) THEN
1107            DO jj=j1,j2
1108               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1109               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
1110               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
1111            END DO
1112         ENDIF
1113         !
1114      ENDIF
1115      !
1116   END SUBROUTINE reflux_sshu
1117
1118   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1119      !!----------------------------------------------------------------------
1120      !!                      *** ROUTINE updatevb2b ***
1121      !!----------------------------------------------------------------------
1122      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1123      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
1124      LOGICAL                         , INTENT(in   ) ::   before
1125      !!
1126      INTEGER :: ji, jj
1127      REAL(wp) :: zrhox, za1, zcor
1128      !!---------------------------------------------
1129      !
1130      IF( before ) THEN
1131         zrhox = Agrif_Rhox()
1132         DO jj=j1,j2
1133            DO ji=i1,i2
1134               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1135            END DO
1136         END DO
1137         tabres = zrhox * tabres
1138      ELSE
1139         !
1140         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1141         !
1142         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1143         DO jj=j1,j2
1144            DO ji=i1,i2
1145               zcor=tabres(ji,jj) - vb2_b(ji,jj)
1146               ! Update time integrated fluxes also in case of multiply nested grids:
1147               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
1148               ! Update corrective fluxes:
1149               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
1150               ! Update half step back fluxes:
1151               vb2_b(ji,jj) = tabres(ji,jj)
1152            END DO
1153         END DO
1154      ENDIF
1155      !
1156   END SUBROUTINE updatevb2b
1157
1158   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1159      !!---------------------------------------------
1160      !!          *** ROUTINE reflux_sshv ***
1161      !!---------------------------------------------
1162      INTEGER, INTENT(in) :: i1, i2, j1, j2
1163      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1164      LOGICAL, INTENT(in) :: before
1165      INTEGER, INTENT(in) :: nb, ndir
1166      !!
1167      LOGICAL :: southern_side, northern_side 
1168      INTEGER :: ji, jj
1169      REAL(wp) :: zrhox, za1, zcor
1170      !!---------------------------------------------
1171      !
1172      IF (before) THEN
1173         zrhox = Agrif_Rhox()
1174         DO jj=j1,j2
1175            DO ji=i1,i2
1176               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1177            END DO
1178         END DO
1179         tabres = zrhox * tabres
1180      ELSE
1181         !
1182         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1183         !
1184         southern_side = (nb == 2).AND.(ndir == 1)
1185         northern_side = (nb == 2).AND.(ndir == 2)
1186         !
1187         IF (southern_side) THEN
1188            DO ji=i1,i2
1189               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
1190               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
1191               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
1192            END DO
1193         ENDIF
1194         IF (northern_side) THEN               
1195            DO ji=i1,i2
1196               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
1197               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
1198               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
1199            END DO
1200         ENDIF
1201         !
1202      ENDIF
1203      !
1204   END SUBROUTINE reflux_sshv
1205
1206   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1207      !
1208      ! ====>>>>>>>>>>    currently not used
1209      !
1210      !!----------------------------------------------------------------------
1211      !!                      *** ROUTINE updateT ***
1212      !!----------------------------------------------------------------------
1213      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
1214      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
1215      LOGICAL                                    , INTENT(in   ) ::   before
1216      !!
1217      INTEGER :: ji,jj,jk
1218      REAL(wp) :: ztemp
1219      !!----------------------------------------------------------------------
1220
1221      IF (before) THEN
1222         DO jk=k1,k2
1223            DO jj=j1,j2
1224               DO ji=i1,i2
1225                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1226                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1227                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1228               END DO
1229            END DO
1230         END DO
1231         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
1232         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1233         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1234      ELSE
1235         DO jk=k1,k2
1236            DO jj=j1,j2
1237               DO ji=i1,i2
1238                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1239                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1240                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
1241                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
1242                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1243                     print *,'CORR = ',ztemp-1.
1244                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
1245                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
1246                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1247                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1248                  END IF
1249               END DO
1250            END DO
1251         END DO
1252      ENDIF
1253      !
1254   END SUBROUTINE update_scales
1255
1256
1257   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1258      !!----------------------------------------------------------------------
1259      !!                      *** ROUTINE updateen ***
1260      !!----------------------------------------------------------------------
1261      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1262      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1263      LOGICAL                               , INTENT(in   ) ::   before
1264      !!----------------------------------------------------------------------
1265      !
1266      IF( before ) THEN
1267         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1268      ELSE
1269         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1270      ENDIF
1271      !
1272   END SUBROUTINE updateEN
1273
1274
1275   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1276      !!----------------------------------------------------------------------
1277      !!                      *** ROUTINE updateavt ***
1278      !!----------------------------------------------------------------------
1279      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1280      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1281      LOGICAL                               , INTENT(in   ) ::   before
1282      !!----------------------------------------------------------------------
1283      !
1284      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1285      ELSE                ;   avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1286      ENDIF
1287      !
1288   END SUBROUTINE updateAVT
1289
1290
1291   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1292      !!---------------------------------------------
1293      !!           *** ROUTINE updateavm ***
1294      !!----------------------------------------------------------------------
1295      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1296      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1297      LOGICAL                               , INTENT(in   ) ::   before
1298      !!----------------------------------------------------------------------
1299      !
1300      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1301      ELSE                ;   avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1302      ENDIF
1303      !
1304   END SUBROUTINE updateAVM
1305
1306   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1307      !!---------------------------------------------
1308      !!           *** ROUTINE updatee3t ***
1309      !!---------------------------------------------
1310      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1311      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1312      LOGICAL, INTENT(in) :: before
1313      !
1314      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1315      INTEGER :: ji,jj,jk
1316      REAL(wp) :: zcoef
1317      !!---------------------------------------------
1318      !
1319      IF (.NOT.before) THEN
1320         !
1321         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) 
1322         !
1323         ! Update e3t from ssh (z* case only)
1324         DO jk = 1, jpkm1
1325            DO jj=j1,j2
1326               DO ji=i1,i2
1327                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
1328                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1329               END DO
1330            END DO
1331         END DO
1332         !
1333         ! 1) Updates at BEFORE time step:
1334         ! -------------------------------
1335         !
1336         ! Save "old" scale factor (prior update) for subsequent asselin correction
1337         ! of prognostic variables
1338         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1339
1340         ! One should also save e3t_b, but lacking of workspace...
1341!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
1342
1343         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1344            DO jk = 1, jpkm1
1345               DO jj=j1,j2
1346                  DO ji=i1,i2
1347                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
1348                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1349                  END DO
1350               END DO
1351            END DO
1352            !
1353            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)
1354            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1355            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1356            !
1357            DO jk = 2, jpk
1358               DO jj = j1,j2
1359                  DO ji = i1,i2           
1360                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1361                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
1362                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
1363                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
1364                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1365                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1366                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
1367                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
1368                  END DO
1369               END DO
1370            END DO
1371            !
1372         ENDIF       
1373         !
1374         ! 2) Updates at NOW time step:
1375         ! ----------------------------
1376         !
1377         ! Update vertical scale factor at T-points:
1378         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1379         !
1380         ! Update total depth:
1381         ht_n(i1:i2,j1:j2) = 0._wp
1382         DO jk = 1, jpkm1
1383            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)
1384         END DO
1385         !
1386         ! Update vertical scale factor at W-points and depths:
1387         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)
1388         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1389         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1390         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
1391         !
1392         DO jk = 2, jpk
1393            DO jj = j1,j2
1394               DO ji = i1,i2           
1395               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1396               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) )   &
1397               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1398               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1399               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
1400                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
1401               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
1402               END DO
1403            END DO
1404         END DO
1405         !
1406         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1407            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
1408            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
1409            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1410            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1411         ENDIF
1412         !
1413         DEALLOCATE(ptab)
1414      ENDIF
1415      !
1416   END SUBROUTINE updatee3t
1417
1418#else
1419   !!----------------------------------------------------------------------
1420   !!   Empty module                                          no AGRIF zoom
1421   !!----------------------------------------------------------------------
1422CONTAINS
1423   SUBROUTINE agrif_oce_update_empty
1424      WRITE(*,*)  'agrif_oce_update : You should not have seen this print! error?'
1425   END SUBROUTINE agrif_oce_update_empty
1426#endif
1427
1428   !!======================================================================
1429END MODULE agrif_oce_update
1430
Note: See TracBrowser for help on using the repository browser.