source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/NST/agrif_oce_update.F90 @ 9863

Last change on this file since 9863 was 9863, checked in by gm, 3 years ago

#1911 (ENHANCE-04): simplified implementation of the Euler stepping at nit000

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