source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/NST/agrif_oce_update.F90 @ 13219

Last change on this file since 13219 was 13219, checked in by smasson, 3 months ago

better e3: update with trunk@13218 see #2385

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