source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_update.F90 @ 13565

Last change on this file since 13565 was 13565, checked in by jchanut, 5 months ago

#2222, 1) Added parent bathymetry volume consistency check 2) Added velocity extrapolation in update 3) Corrected bdy issue #2519

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