source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 10425

Last change on this file since 10425 was 10425, checked in by smasson, 23 months ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 77.1 KB
Line 
1MODULE dynspg_ts
2
3   !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
4
5   !!======================================================================
6   !!                   ***  MODULE  dynspg_ts  ***
7   !! Ocean dynamics:  surface pressure gradient trend, split-explicit scheme
8   !!======================================================================
9   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
10   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
11   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
12   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
13   !!              -   ! 2008-01  (R. Benshila)  change averaging method
14   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
15   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
17   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
18   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
19   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
20   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
21   !!             4.0  ! 2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
22   !!---------------------------------------------------------------------
23
24   !!----------------------------------------------------------------------
25   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
26   !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
27   !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
28   !!   ts_rst         : read/write time-splitting fields in restart file
29   !!----------------------------------------------------------------------
30   USE oce             ! ocean dynamics and tracers
31   USE dom_oce         ! ocean space and time domain
32   USE sbc_oce         ! surface boundary condition: ocean
33   USE zdf_oce         ! vertical physics: variables
34   USE zdfdrg          ! vertical physics: top/bottom drag coef.
35   USE sbcisf          ! ice shelf variable (fwfisf)
36   USE sbcapr          ! surface boundary condition: atmospheric pressure
37   USE dynadv    , ONLY: ln_dynadv_vec
38   USE dynvor          ! vortivity scheme indicators
39   USE phycst          ! physical constants
40   USE dynvor          ! vorticity term
41   USE wet_dry         ! wetting/drying flux limter
42   USE bdy_oce         ! open boundary
43   USE bdytides        ! open boundary condition data
44   USE bdydyn2d        ! open boundary conditions on barotropic variables
45   USE sbctide         ! tides
46   USE updtide         ! tide potential
47   USE sbcwave         ! surface wave
48   USE diatmb          ! Top,middle,bottom output
49#if defined key_agrif
50   USE agrif_oce_interp ! agrif
51   USE agrif_oce
52#endif
53#if defined key_asminc   
54   USE asminc          ! Assimilation increment
55#endif
56   !
57   USE in_out_manager  ! I/O manager
58   USE lib_mpp         ! distributed memory computing library
59   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl          ! Print control
61   USE iom             ! IOM library
62   USE restart         ! only for lrst_oce
63   USE diatmb          ! Top,middle,bottom output
64
65   IMPLICIT NONE
66   PRIVATE
67
68   PUBLIC dyn_spg_ts        ! called by dyn_spg
69   PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
70
71   !! Time filtered arrays at baroclinic time step:
72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step
73   !
74   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
75   REAL(wp),SAVE :: rdtbt       ! Barotropic time step
76   !
77   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields
78   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwz                ! ff_f/h at F points
79   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftnw, ftne         ! triad of coriolis parameter
80   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftsw, ftse         ! (only used with een vorticity scheme)
81
82   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios
83   REAL(wp) ::   r1_8  = 0.125_wp         !
84   REAL(wp) ::   r1_4  = 0.25_wp          !
85   REAL(wp) ::   r1_2  = 0.5_wp           !
86
87   !! * Substitutions
88#  include "vectopt_loop_substitute.h90"
89   !!----------------------------------------------------------------------
90   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
91   !! $Id$
92   !! Software governed by the CeCILL license (see ./LICENSE)
93   !!----------------------------------------------------------------------
94CONTAINS
95
96   INTEGER FUNCTION dyn_spg_ts_alloc()
97      !!----------------------------------------------------------------------
98      !!                  ***  routine dyn_spg_ts_alloc  ***
99      !!----------------------------------------------------------------------
100      INTEGER :: ierr(3)
101      !!----------------------------------------------------------------------
102      ierr(:) = 0
103      !
104      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
105      !
106      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
107         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
108         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
109         !
110      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
111      !
112      dyn_spg_ts_alloc = MAXVAL( ierr(:) )
113      !
114      CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
115      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
116      !
117   END FUNCTION dyn_spg_ts_alloc
118
119
120   SUBROUTINE dyn_spg_ts( kt )
121      !!----------------------------------------------------------------------
122      !!
123      !! ** Purpose : - Compute the now trend due to the explicit time stepping
124      !!              of the quasi-linear barotropic system, and add it to the
125      !!              general momentum trend.
126      !!
127      !! ** Method  : - split-explicit schem (time splitting) :
128      !!      Barotropic variables are advanced from internal time steps
129      !!      "n"   to "n+1" if ln_bt_fw=T
130      !!      or from
131      !!      "n-1" to "n+1" if ln_bt_fw=F
132      !!      thanks to a generalized forward-backward time stepping (see ref. below).
133      !!
134      !! ** Action :
135      !!      -Update the filtered free surface at step "n+1"      : ssha
136      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
137      !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
138      !!      These are used to advect tracers and are compliant with discrete
139      !!      continuity equation taken at the baroclinic time steps. This
140      !!      ensures tracers conservation.
141      !!      - (ua, va) momentum trend updated with barotropic component.
142      !!
143      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
144      !!---------------------------------------------------------------------
145      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
146      !
147      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
148      LOGICAL  ::   ll_fw_start           ! =T : forward integration
149      LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
150      LOGICAL  ::   ll_tmp1, ll_tmp2      ! local logical variables used in W/D
151      INTEGER  ::   ikbu, iktu, noffset   ! local integers
152      INTEGER  ::   ikbv, iktv            !   -      -
153      REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars
154      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      -
155      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      -
156      REAL(wp) ::   za0, za1, za2, za3              !   -      -
157      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      -
158      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
159      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
160      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
161      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e
162      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
163      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
164      !
165      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.
166
167      INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point
168
169      REAL(wp) ::   zepsilon, zgamma            !   -      -
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy   ! Wetting/Dying gravity filter coef.
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
173      !!----------------------------------------------------------------------
174      !
175      IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
176      !                                         !* Allocate temporary arrays
177      IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
178      !
179      zmdi=1.e+20                               !  missing data indicator for masking
180      !
181      zwdramp = r_rn_wdmin1               ! simplest ramp
182!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
183      !                                         ! reciprocal of baroclinic time step
184      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
185      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt
186      ENDIF
187      r1_2dt_b = 1.0_wp / z2dt_bf 
188      !
189      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart
190      ll_fw_start = .FALSE.
191      !                                         ! time offset in steps for bdy data update
192      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro
193      ELSE                       ;   noffset =   0 
194      ENDIF
195      !
196      IF( kt == nit000 ) THEN                   !* initialisation
197         !
198         IF(lwp) WRITE(numout,*)
199         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
200         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
201         IF(lwp) WRITE(numout,*)
202         !
203         IF( neuler == 0 )   ll_init=.TRUE.
204         !
205         IF( ln_bt_fw .OR. neuler == 0 ) THEN
206            ll_fw_start =.TRUE.
207            noffset     = 0
208         ELSE
209            ll_fw_start =.FALSE.
210         ENDIF
211         !
212         ! Set averaging weights and cycle length:
213         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
214         !
215      ENDIF
216      !
217      IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities)
218         DO jj = 2, jpjm1
219            DO ji = fs_2, fs_jpim1   ! vector opt.
220               zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
221               zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
222            END DO
223         END DO
224      ELSE                    ! bottom friction only
225         DO jj = 2, jpjm1
226            DO ji = fs_2, fs_jpim1   ! vector opt.
227               zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
228               zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
229            END DO
230         END DO
231      ENDIF
232      !
233      ! Set arrays to remove/compute coriolis trend.
234      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
235      ! Note that these arrays are also used during barotropic loop. These are however frozen
236      ! although they should be updated in the variable volume case. Not a big approximation.
237      ! To remove this approximation, copy lines below inside barotropic loop
238      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
239      !
240      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
241         !
242         SELECT CASE( nvor_scheme )
243         CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
244            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
245            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
246               DO jj = 1, jpjm1
247                  DO ji = 1, jpim1
248                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
249                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp 
250                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
251                  END DO
252               END DO
253            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
254               DO jj = 1, jpjm1
255                  DO ji = 1, jpim1
256                     zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
257                        &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
258                        &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
259                        &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
260                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
261                  END DO
262               END DO
263            END SELECT
264            CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
265            !
266            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
267            DO jj = 2, jpj
268               DO ji = 2, jpi
269                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
270                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
271                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
272                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
273               END DO
274            END DO
275            !
276         CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
277            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
278            DO jj = 2, jpj
279               DO ji = 2, jpi
280                  z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
281                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
282                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
283                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
284                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
285               END DO
286            END DO
287            !
288         CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
289            !
290            zwz(:,:) = 0._wp
291            zhf(:,:) = 0._wp
292           
293!!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
294!!gm    A priori a better value should be something like :
295!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
296!!gm                     divided by the sum of the corresponding mask
297!!gm
298!!           
299            IF( .NOT.ln_sco ) THEN
300 
301   !!gm  agree the JC comment  : this should be done in a much clear way
302 
303   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
304   !     Set it to zero for the time being
305   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
306   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
307   !              ENDIF
308   !              zhf(:,:) = gdepw_0(:,:,jk+1)
309               !
310            ELSE
311               !
312               !zhf(:,:) = hbatf(:,:)
313               DO jj = 1, jpjm1
314                  DO ji = 1, jpim1
315                     zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
316                        &              + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
317                        &       / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
318                        &              + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
319                  END DO
320               END DO
321            ENDIF
322            !
323            DO jj = 1, jpjm1
324               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
325            END DO
326            !
327            DO jk = 1, jpkm1
328               DO jj = 1, jpjm1
329                  zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
330               END DO
331            END DO
332            CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
333            ! JC: TBC. hf should be greater than 0
334            DO jj = 1, jpj
335               DO ji = 1, jpi
336                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
337               END DO
338            END DO
339            zwz(:,:) = ff_f(:,:) * zwz(:,:)
340         END SELECT
341      ENDIF
342      !
343      ! If forward start at previous time step, and centered integration,
344      ! then update averaging weights:
345      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
346         ll_fw_start=.FALSE.
347         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
348      ENDIF
349                         
350      ! -----------------------------------------------------------------------------
351      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
352      ! -----------------------------------------------------------------------------
353      !     
354      !
355      !                                   !* e3*d/dt(Ua) (Vertically integrated)
356      !                                   ! --------------------------------------------------
357      zu_frc(:,:) = 0._wp
358      zv_frc(:,:) = 0._wp
359      !
360      DO jk = 1, jpkm1
361         zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
362         zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
363      END DO
364      !
365      zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
366      zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
367      !
368      !
369      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
370      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
371         DO jj = 2, jpjm1
372            DO ji = fs_2, fs_jpim1   ! vector opt.
373               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
374               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
375            END DO
376         END DO
377      END DO
378     
379!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
380!!gm  Is it correct to do so ?   I think so...
381     
382     
383      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
384      !                                   ! --------------------------------------------------------
385      !
386      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes
387      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
388      !
389      SELECT CASE( nvor_scheme )
390      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
391         DO jj = 2, jpjm1
392            DO ji = 2, jpim1   ! vector opt.
393               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    &
394                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
395                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
396                  !
397               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    &
398                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
399                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
400            END DO 
401         END DO 
402         !         
403      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
404         DO jj = 2, jpjm1
405            DO ji = fs_2, fs_jpim1   ! vector opt.
406               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
407               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
408               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
409               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
410               ! energy conserving formulation for planetary vorticity term
411               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
412               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
413            END DO
414         END DO
415         !
416      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
417         DO jj = 2, jpjm1
418            DO ji = fs_2, fs_jpim1   ! vector opt.
419               zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
420                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
421               zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
422                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
423               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
424               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
425            END DO
426         END DO
427         !
428      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
429         DO jj = 2, jpjm1
430            DO ji = fs_2, fs_jpim1   ! vector opt.
431               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
432                &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
433                &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
434                &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
435               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
436                &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
437                &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
438                &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
439            END DO
440         END DO
441         !
442      END SELECT
443      !
444      !                                   !* Right-Hand-Side of the barotropic momentum equation
445      !                                   ! ----------------------------------------------------
446      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
447         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters
448            DO jj = 2, jpjm1
449               DO ji = 2, jpim1 
450                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
451                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
452                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
453                     &                                                         > rn_wdmin1 + rn_wdmin2
454                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
455                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
456                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
457                  IF(ll_tmp1) THEN
458                     zcpx(ji,jj) = 1.0_wp
459                  ELSEIF(ll_tmp2) THEN
460                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
461                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
462                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
463                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
464                  ELSE
465                     zcpx(ji,jj) = 0._wp
466                  ENDIF
467                  !
468                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
469                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
470                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
471                     &                                                       > rn_wdmin1 + rn_wdmin2
472                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
473                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
474                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
475 
476                  IF(ll_tmp1) THEN
477                     zcpy(ji,jj) = 1.0_wp
478                  ELSE IF(ll_tmp2) THEN
479                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
480                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
481                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
482                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
483                  ELSE
484                     zcpy(ji,jj) = 0._wp
485                  ENDIF
486               END DO
487            END DO
488            !
489            DO jj = 2, jpjm1
490               DO ji = 2, jpim1
491                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
492                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
493                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
494                     &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
495               END DO
496            END DO
497            !
498         ELSE
499            !
500            DO jj = 2, jpjm1
501               DO ji = fs_2, fs_jpim1   ! vector opt.
502                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
503                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
504               END DO
505            END DO
506         ENDIF
507         !
508      ENDIF
509      !
510      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
511         DO ji = fs_2, fs_jpim1
512             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
513             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
514          END DO
515      END DO 
516      !
517      !                                         ! Add bottom stress contribution from baroclinic velocities:     
518      IF (ln_bt_fw) THEN
519         DO jj = 2, jpjm1                         
520            DO ji = fs_2, fs_jpim1   ! vector opt.
521               ikbu = mbku(ji,jj)       
522               ikbv = mbkv(ji,jj)   
523               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
524               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
525            END DO
526         END DO
527      ELSE
528         DO jj = 2, jpjm1
529            DO ji = fs_2, fs_jpim1   ! vector opt.
530               ikbu = mbku(ji,jj)       
531               ikbv = mbkv(ji,jj)   
532               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
533               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
534            END DO
535         END DO
536      ENDIF
537      !
538      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
539      IF( ln_wd_il ) THEN
540         zztmp = -1._wp / rdtbt
541         DO jj = 2, jpjm1
542            DO ji = fs_2, fs_jpim1   ! vector opt.
543               zu_frc(ji,jj) = zu_frc(ji,jj) + & 
544               & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) *  wdrampu(ji,jj)
545               zv_frc(ji,jj) = zv_frc(ji,jj) + & 
546               & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) *  wdrampv(ji,jj)
547            END DO
548         END DO
549      ELSE
550         DO jj = 2, jpjm1
551            DO ji = fs_2, fs_jpim1   ! vector opt.
552               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)
553               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)
554            END DO
555         END DO
556      END IF
557      !
558      IF( ln_isfcav ) THEN       ! Add TOP stress contribution from baroclinic velocities:     
559         IF( ln_bt_fw ) THEN
560            DO jj = 2, jpjm1
561               DO ji = fs_2, fs_jpim1   ! vector opt.
562                  iktu = miku(ji,jj)
563                  iktv = mikv(ji,jj)
564                  zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
565                  zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
566               END DO
567            END DO
568         ELSE
569            DO jj = 2, jpjm1
570               DO ji = fs_2, fs_jpim1   ! vector opt.
571                  iktu = miku(ji,jj)
572                  iktv = mikv(ji,jj)
573                  zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
574                  zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
575               END DO
576            END DO
577         ENDIF
578         !
579         ! Note that the "unclipped" top friction parameter is used even with explicit drag
580         DO jj = 2, jpjm1             
581            DO ji = fs_2, fs_jpim1   ! vector opt.
582               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)
583               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)
584            END DO
585         END DO
586      ENDIF
587      !       
588      IF( ln_bt_fw ) THEN                        ! Add wind forcing
589         DO jj = 2, jpjm1
590            DO ji = fs_2, fs_jpim1   ! vector opt.
591               zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
592               zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
593            END DO
594         END DO
595      ELSE
596         zztmp = r1_rau0 * r1_2
597         DO jj = 2, jpjm1
598            DO ji = fs_2, fs_jpim1   ! vector opt.
599               zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
600               zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
601            END DO
602         END DO
603      ENDIF 
604      !
605      IF( ln_apr_dyn ) THEN                     ! Add atm pressure forcing
606         IF( ln_bt_fw ) THEN
607            DO jj = 2, jpjm1             
608               DO ji = fs_2, fs_jpim1   ! vector opt.
609                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
610                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
611                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
612                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
613               END DO
614            END DO
615         ELSE
616            zztmp = grav * r1_2
617            DO jj = 2, jpjm1             
618               DO ji = fs_2, fs_jpim1   ! vector opt.
619                  zu_spg = zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
620                      &             + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
621                  zv_spg = zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
622                      &             + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
623                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
624                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
625               END DO
626            END DO
627         ENDIF
628      ENDIF
629      !                                   !* Right-Hand-Side of the barotropic ssh equation
630      !                                   ! -----------------------------------------------
631      !                                         ! Surface net water flux and rivers
632      IF (ln_bt_fw) THEN
633         zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
634      ELSE
635         zztmp = r1_rau0 * r1_2
636         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
637                &                 + fwfisf(:,:) + fwfisf_b(:,:)                     )
638      ENDIF
639      !
640      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary
641         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
642      ENDIF
643      !
644#if defined key_asminc
645      !                                         ! Include the IAU weighted SSH increment
646      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
647         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
648      ENDIF
649#endif
650      !                                   !* Fill boundary data arrays for AGRIF
651      !                                   ! ------------------------------------
652#if defined key_agrif
653         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
654#endif
655      !
656      ! -----------------------------------------------------------------------
657      !  Phase 2 : Integration of the barotropic equations
658      ! -----------------------------------------------------------------------
659      !
660      !                                             ! ==================== !
661      !                                             !    Initialisations   !
662      !                                             ! ==================== ! 
663      ! Initialize barotropic variables:     
664      IF( ll_init )THEN
665         sshbb_e(:,:) = 0._wp
666         ubb_e  (:,:) = 0._wp
667         vbb_e  (:,:) = 0._wp
668         sshb_e (:,:) = 0._wp
669         ub_e   (:,:) = 0._wp
670         vb_e   (:,:) = 0._wp
671      ENDIF
672
673      !
674      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
675         sshn_e(:,:) =    sshn(:,:)           
676         un_e  (:,:) =    un_b(:,:)           
677         vn_e  (:,:) =    vn_b(:,:)
678         !
679         hu_e  (:,:) =    hu_n(:,:)       
680         hv_e  (:,:) =    hv_n(:,:) 
681         hur_e (:,:) = r1_hu_n(:,:)   
682         hvr_e (:,:) = r1_hv_n(:,:)
683      ELSE                                ! CENTRED integration: start from BEFORE fields
684         sshn_e(:,:) =    sshb(:,:)
685         un_e  (:,:) =    ub_b(:,:)         
686         vn_e  (:,:) =    vb_b(:,:)
687         !
688         hu_e  (:,:) =    hu_b(:,:)       
689         hv_e  (:,:) =    hv_b(:,:) 
690         hur_e (:,:) = r1_hu_b(:,:)   
691         hvr_e (:,:) = r1_hv_b(:,:)
692      ENDIF
693      !
694      !
695      !
696      ! Initialize sums:
697      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
698      va_b  (:,:) = 0._wp
699      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
700      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
701      vn_adv(:,:) = 0._wp
702      !
703      IF( ln_wd_dl ) THEN
704         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)
705         zvwdmask(:,:) = 0._wp  !
706         zuwdav2 (:,:) = 0._wp 
707         zvwdav2 (:,:) = 0._wp   
708      END IF 
709
710      !                                             ! ==================== !
711      DO jn = 1, icycle                             !  sub-time-step loop  !
712         !                                          ! ==================== !
713         !
714         l_full_nf_update = jn == icycle   ! false: disable full North fold update (performances) for jn = 1 to icycle-1
715         !                                                !  ------------------
716         !                                                !* Update the forcing (BDY and tides)
717         !                                                !  ------------------
718         ! Update only tidal forcing at open boundaries
719         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
720         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
721         !
722         ! Set extrapolation coefficients for predictor step:
723         IF ((jn<3).AND.ll_init) THEN      ! Forward           
724           za1 = 1._wp                                         
725           za2 = 0._wp                       
726           za3 = 0._wp                       
727         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
728           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
729           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
730           za3 =  0.281105_wp              ! za3 = bet
731         ENDIF
732
733         ! Extrapolate barotropic velocities at step jit+0.5:
734         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
735         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
736
737         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
738            !                                             !  ------------------
739            ! Extrapolate Sea Level at step jit+0.5:
740            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
741           
742            ! set wetting & drying mask at tracer points for this barotropic sub-step
743            IF ( ln_wd_dl ) THEN 
744               !
745               IF ( ln_wd_dl_rmp ) THEN
746                  DO jj = 1, jpj                                 
747                     DO ji = 1, jpi   ! vector opt. 
748                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
749!                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN
750                           ztwdmask(ji,jj) = 1._wp
751                        ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
752                           ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) ) 
753                        ELSE
754                           ztwdmask(ji,jj) = 0._wp
755                        END IF
756                     END DO
757                  END DO
758               ELSE
759                  DO jj = 1, jpj                                 
760                     DO ji = 1, jpi   ! vector opt. 
761                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
762                           ztwdmask(ji,jj) = 1._wp
763                        ELSE
764                           ztwdmask(ji,jj) = 0._wp
765                        ENDIF
766                     END DO
767                  END DO
768               ENDIF 
769               !
770            ENDIF 
771            !
772            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
773               DO ji = 2, fs_jpim1   ! Vector opt.
774                  zwx(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
775                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
776                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
777                  zwy(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
778                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
779                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
780               END DO
781            END DO
782            CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp )
783            !
784            zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
785            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:)
786            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
787         ELSE
788            zhup2_e(:,:) = hu_n(:,:)
789            zhvp2_e(:,:) = hv_n(:,:)
790            zhtp2_e(:,:) = ht_n(:,:)
791         ENDIF
792         !                                                !* after ssh
793         !                                                !  -----------
794         ! One should enforce volume conservation at open boundaries here
795         ! considering fluxes below:
796         !
797         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
798         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
799         !
800#if defined key_agrif
801         ! Set fluxes during predictor step to ensure volume conservation
802         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
803            IF((nbondi == -1).OR.(nbondi == 2)) THEN
804               DO jj = 1, jpj
805                  zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
806                  zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
807               END DO
808            ENDIF
809            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
810               DO jj=1,jpj
811                  zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
812                  zwy(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj)
813               END DO
814            ENDIF
815            IF((nbondj == -1).OR.(nbondj == 2)) THEN
816               DO ji=1,jpi
817                  zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
818                  zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
819               END DO
820            ENDIF
821            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
822               DO ji=1,jpi
823                  zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
824                  zwx(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1)
825               END DO
826            ENDIF
827         ENDIF
828#endif
829         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
830
831         IF ( ln_wd_dl ) THEN 
832            !
833            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
834            !
835            DO jj = 1, jpjm1                                 
836               DO ji = 1, jpim1   
837                  IF ( zwx(ji,jj) > 0.0 ) THEN
838                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj) 
839                  ELSE
840                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 
841                  END IF
842                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
843                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
844
845                  IF ( zwy(ji,jj) > 0.0 ) THEN
846                     zvwdmask(ji, jj) = ztwdmask(ji, jj  )
847                  ELSE
848                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 
849                  END IF
850                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 
851                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
852               END DO
853            END DO
854            !
855         ENDIF   
856         
857         ! Sum over sub-time-steps to compute advective velocities
858         za2 = wgtbtp2(jn)
859         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
860         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
861         
862         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
863         IF ( ln_wd_dl_bc ) THEN
864            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
865            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
866         END IF 
867
868         ! Set next sea level:
869         DO jj = 2, jpjm1                                 
870            DO ji = fs_2, fs_jpim1   ! vector opt.
871               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
872                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
873            END DO
874         END DO
875         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
876         
877         CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T',  1._wp )
878
879         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
880         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
881#if defined key_agrif
882         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
883#endif
884         
885         ! Sea Surface Height at u-,v-points (vvl case only)
886         IF( .NOT.ln_linssh ) THEN                               
887            DO jj = 2, jpjm1
888               DO ji = 2, jpim1      ! NO Vector Opt.
889                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
890                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
891                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
892                  zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
893                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
894                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
895               END DO
896            END DO
897            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
898         ENDIF   
899         !                                 
900         ! Half-step back interpolation of SSH for surface pressure computation:
901         !----------------------------------------------------------------------
902         IF ((jn==1).AND.ll_init) THEN
903           za0=1._wp                        ! Forward-backward
904           za1=0._wp                           
905           za2=0._wp
906           za3=0._wp
907         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
908           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
909           za1=-0.1666666666666_wp          ! za1 = gam
910           za2= 0.0833333333333_wp          ! za2 = eps
911           za3= 0._wp             
912         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
913            IF (rn_bt_alpha==0._wp) THEN
914               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps
915               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps
916               za2=0.088_wp                 ! za2 = gam
917               za3=0.013_wp                 ! za3 = eps
918            ELSE
919               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
920               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
921               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
922               za1 = 1._wp - za0 - zgamma - zepsilon
923               za2 = zgamma
924               za3 = zepsilon
925            ENDIF
926         ENDIF
927         !
928         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
929            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
930         
931         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters
932           DO jj = 2, jpjm1
933              DO ji = 2, jpim1 
934                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
935                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            &
936                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
937                     &                                                             > rn_wdmin1 + rn_wdmin2
938                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
939                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
940                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
941   
942                IF(ll_tmp1) THEN
943                  zcpx(ji,jj) = 1.0_wp
944                ELSE IF(ll_tmp2) THEN
945                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here
946                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
947                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) )
948                ELSE
949                  zcpx(ji,jj) = 0._wp
950                ENDIF
951                !
952                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
953                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            &
954                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
955                     &                                                             > rn_wdmin1 + rn_wdmin2
956                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
957                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
958                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
959   
960                IF(ll_tmp1) THEN
961                  zcpy(ji,jj) = 1.0_wp
962                ELSEIF(ll_tmp2) THEN
963                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
964                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
965                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) )
966                ELSE
967                  zcpy(ji,jj) = 0._wp
968                ENDIF
969              END DO
970           END DO
971         ENDIF
972         !
973         ! Compute associated depths at U and V points:
974         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form
975            !                                       
976            DO jj = 2, jpjm1                           
977               DO ji = 2, jpim1
978                  zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
979                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
980                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
981                  zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
982                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
983                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
984                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
985                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
986               END DO
987            END DO
988            !
989         ENDIF
990         !
991         ! Add Coriolis trend:
992         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
993         ! at each time step. We however keep them constant here for optimization.
994         ! Recall that zwx and zwy arrays hold fluxes at this stage:
995         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
996         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
997         !
998         SELECT CASE( nvor_scheme )
999         CASE( np_ENT )             ! energy conserving scheme (t-point)
1000         DO jj = 2, jpjm1
1001            DO ji = 2, jpim1   ! vector opt.
1002
1003               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
1004               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
1005           
1006               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   &
1007                  &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   &
1008                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   )
1009                  !
1010               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1011                  &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   & 
1012                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   ) 
1013            END DO 
1014         END DO 
1015         !         
1016         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point)
1017            DO jj = 2, jpjm1
1018               DO ji = fs_2, fs_jpim1   ! vector opt.
1019                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
1020                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
1021                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
1022                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
1023                  zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1024                  zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1025               END DO
1026            END DO
1027            !
1028         CASE( np_ENS )             ! enstrophy conserving scheme (f-point)
1029            DO jj = 2, jpjm1
1030               DO ji = fs_2, fs_jpim1   ! vector opt.
1031                  zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
1032                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
1033                  zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
1034                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
1035                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1036                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1037               END DO
1038            END DO
1039            !
1040         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f)
1041            DO jj = 2, jpjm1
1042               DO ji = fs_2, fs_jpim1   ! vector opt.
1043                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  &
1044                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  &
1045                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  & 
1046                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  )
1047                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  & 
1048                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  &
1049                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  & 
1050                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  )
1051               END DO
1052            END DO
1053            !
1054         END SELECT
1055         !
1056         ! Add tidal astronomical forcing if defined
1057         IF ( ln_tide .AND. ln_tide_pot ) THEN
1058            DO jj = 2, jpjm1
1059               DO ji = fs_2, fs_jpim1   ! vector opt.
1060                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1061                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1062                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1063                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1064               END DO
1065            END DO
1066         ENDIF
1067         !
1068         ! Add bottom stresses:
1069!jth do implicitly instead
1070         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1071            DO jj = 2, jpjm1
1072               DO ji = fs_2, fs_jpim1   ! vector opt.
1073                  zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1074                  zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1075               END DO
1076            END DO
1077         ENDIF 
1078         !
1079         ! Surface pressure trend:
1080         IF( ln_wd_il ) THEN
1081           DO jj = 2, jpjm1
1082              DO ji = 2, jpim1 
1083                 ! Add surface pressure gradient
1084                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1085                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1086                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 
1087                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1088              END DO
1089           END DO
1090         ELSE
1091           DO jj = 2, jpjm1
1092              DO ji = fs_2, fs_jpim1   ! vector opt.
1093                 ! Add surface pressure gradient
1094                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1095                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1096                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1097                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1098              END DO
1099           END DO
1100         END IF
1101
1102         !
1103         ! Set next velocities:
1104         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
1105            DO jj = 2, jpjm1
1106               DO ji = fs_2, fs_jpim1   ! vector opt.
1107                  ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
1108                            &     + rdtbt * (                      zwx(ji,jj)   &
1109                            &                                 + zu_trd(ji,jj)   &
1110                            &                                 + zu_frc(ji,jj) ) & 
1111                            &   ) * ssumask(ji,jj)
1112
1113                  va_e(ji,jj) = (                                 vn_e(ji,jj)   &
1114                            &     + rdtbt * (                      zwy(ji,jj)   &
1115                            &                                 + zv_trd(ji,jj)   &
1116                            &                                 + zv_frc(ji,jj) ) &
1117                            &   ) * ssvmask(ji,jj)
1118 
1119               END DO
1120            END DO
1121            !
1122         ELSE                           !* Flux form
1123            DO jj = 2, jpjm1
1124               DO ji = fs_2, fs_jpim1   ! vector opt.
1125
1126                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1127                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1128
1129                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1130                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1131
1132                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
1133                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
1134                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
1135                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
1136                            &   ) * zhura
1137
1138                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
1139                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
1140                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
1141                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
1142                            &   ) * zhvra
1143               END DO
1144            END DO
1145         ENDIF
1146!jth implicit bottom friction:
1147         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1148            DO jj = 2, jpjm1
1149               DO ji = fs_2, fs_jpim1   ! vector opt.
1150                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1151                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1152               END DO
1153            END DO
1154         ENDIF
1155
1156         
1157         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
1158            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1159            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1160            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1161            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1162            !
1163         ENDIF
1164         !                                             !* domain lateral boundary
1165         CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1166         !
1167         !                                                 ! open boundaries
1168         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1169#if defined key_agrif                                                           
1170         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
1171#endif
1172         !                                             !* Swap
1173         !                                             !  ----
1174         ubb_e  (:,:) = ub_e  (:,:)
1175         ub_e   (:,:) = un_e  (:,:)
1176         un_e   (:,:) = ua_e  (:,:)
1177         !
1178         vbb_e  (:,:) = vb_e  (:,:)
1179         vb_e   (:,:) = vn_e  (:,:)
1180         vn_e   (:,:) = va_e  (:,:)
1181         !
1182         sshbb_e(:,:) = sshb_e(:,:)
1183         sshb_e (:,:) = sshn_e(:,:)
1184         sshn_e (:,:) = ssha_e(:,:)
1185
1186         !                                             !* Sum over whole bt loop
1187         !                                             !  ----------------------
1188         za1 = wgtbtp1(jn)                                   
1189         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
1190            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
1191            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
1192         ELSE                                       ! Sum transports
1193            IF ( .NOT.ln_wd_dl ) THEN 
1194               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
1195               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
1196            ELSE
1197               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
1198               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
1199            END IF
1200         ENDIF
1201         !                                          ! Sum sea level
1202         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1203
1204         !                                                 ! ==================== !
1205      END DO                                               !        end loop      !
1206      !                                                    ! ==================== !
1207      ! -----------------------------------------------------------------------------
1208      ! Phase 3. update the general trend with the barotropic trend
1209      ! -----------------------------------------------------------------------------
1210      !
1211      ! Set advection velocity correction:
1212      IF (ln_bt_fw) THEN
1213         zwx(:,:) = un_adv(:,:)
1214         zwy(:,:) = vn_adv(:,:)
1215         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1216            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1217            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1218            !
1219            ! Update corrective fluxes for next time step:
1220            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1221            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1222         ELSE
1223            un_bf(:,:) = 0._wp
1224            vn_bf(:,:) = 0._wp 
1225         END IF         
1226         ! Save integrated transport for next computation
1227         ub2_b(:,:) = zwx(:,:)
1228         vb2_b(:,:) = zwy(:,:)
1229      ENDIF
1230
1231
1232      !
1233      ! Update barotropic trend:
1234      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1235         DO jk=1,jpkm1
1236            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1237            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1238         END DO
1239      ELSE
1240         ! At this stage, ssha has been corrected: compute new depths at velocity points
1241         DO jj = 1, jpjm1
1242            DO ji = 1, jpim1      ! NO Vector Opt.
1243               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
1244                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      &
1245                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1246               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
1247                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
1248                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1249            END DO
1250         END DO
1251         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1252         !
1253         DO jk=1,jpkm1
1254            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1255            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1256         END DO
1257         ! Save barotropic velocities not transport:
1258         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1259         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1260      ENDIF
1261
1262
1263      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
1264      DO jk = 1, jpkm1
1265         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)
1266         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1267      END DO
1268
1269      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1270         DO jk = 1, jpkm1
1271            un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1272                       & + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk) 
1273            vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) & 
1274                       & + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk) 
1275         END DO
1276      END IF
1277
1278     
1279      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current
1280      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current
1281      !
1282#if defined key_agrif
1283      ! Save time integrated fluxes during child grid integration
1284      ! (used to update coarse grid transports at next time step)
1285      !
1286      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1287         IF( Agrif_NbStepint() == 0 ) THEN
1288            ub2_i_b(:,:) = 0._wp
1289            vb2_i_b(:,:) = 0._wp
1290         END IF
1291         !
1292         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1293         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1294         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1295      ENDIF
1296#endif     
1297      !                                   !* write time-spliting arrays in the restart
1298      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
1299      !
1300      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
1301      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1302      !
1303      IF( ln_diatmb ) THEN
1304         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
1305         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
1306      ENDIF
1307      !
1308   END SUBROUTINE dyn_spg_ts
1309
1310
1311   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1312      !!---------------------------------------------------------------------
1313      !!                   ***  ROUTINE ts_wgt  ***
1314      !!
1315      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1316      !!----------------------------------------------------------------------
1317      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
1318      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
1319      INTEGER, INTENT(inout) :: jpit      ! cycle length   
1320      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
1321                                                         zwgt2    ! Secondary weights
1322     
1323      INTEGER ::  jic, jn, ji                      ! temporary integers
1324      REAL(wp) :: za1, za2
1325      !!----------------------------------------------------------------------
1326
1327      zwgt1(:) = 0._wp
1328      zwgt2(:) = 0._wp
1329
1330      ! Set time index when averaged value is requested
1331      IF (ll_fw) THEN
1332         jic = nn_baro
1333      ELSE
1334         jic = 2 * nn_baro
1335      ENDIF
1336
1337      ! Set primary weights:
1338      IF (ll_av) THEN
1339           ! Define simple boxcar window for primary weights
1340           ! (width = nn_baro, centered around jic)     
1341         SELECT CASE ( nn_bt_flt )
1342              CASE( 0 )  ! No averaging
1343                 zwgt1(jic) = 1._wp
1344                 jpit = jic
1345
1346              CASE( 1 )  ! Boxcar, width = nn_baro
1347                 DO jn = 1, 3*nn_baro
1348                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1349                    IF (za1 < 0.5_wp) THEN
1350                      zwgt1(jn) = 1._wp
1351                      jpit = jn
1352                    ENDIF
1353                 ENDDO
1354
1355              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
1356                 DO jn = 1, 3*nn_baro
1357                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1358                    IF (za1 < 1._wp) THEN
1359                      zwgt1(jn) = 1._wp
1360                      jpit = jn
1361                    ENDIF
1362                 ENDDO
1363              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1364         END SELECT
1365
1366      ELSE ! No time averaging
1367         zwgt1(jic) = 1._wp
1368         jpit = jic
1369      ENDIF
1370   
1371      ! Set secondary weights
1372      DO jn = 1, jpit
1373        DO ji = jn, jpit
1374             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1375        END DO
1376      END DO
1377
1378      ! Normalize weigths:
1379      za1 = 1._wp / SUM(zwgt1(1:jpit))
1380      za2 = 1._wp / SUM(zwgt2(1:jpit))
1381      DO jn = 1, jpit
1382        zwgt1(jn) = zwgt1(jn) * za1
1383        zwgt2(jn) = zwgt2(jn) * za2
1384      END DO
1385      !
1386   END SUBROUTINE ts_wgt
1387
1388
1389   SUBROUTINE ts_rst( kt, cdrw )
1390      !!---------------------------------------------------------------------
1391      !!                   ***  ROUTINE ts_rst  ***
1392      !!
1393      !! ** Purpose : Read or write time-splitting arrays in restart file
1394      !!----------------------------------------------------------------------
1395      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1396      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1397      !!----------------------------------------------------------------------
1398      !
1399      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1400         !                                   ! ---------------
1401         IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN    !* Read the restart file
1402            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
1403            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
1404            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
1405            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
1406            IF( .NOT.ln_bt_av ) THEN
1407               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
1408               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
1409               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
1410               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
1411               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
1412               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
1413            ENDIF
1414#if defined key_agrif
1415            ! Read time integrated fluxes
1416            IF ( .NOT.Agrif_Root() ) THEN
1417               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
1418               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
1419            ENDIF
1420#endif
1421         ELSE                                   !* Start from rest
1422            IF(lwp) WRITE(numout,*)
1423            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
1424            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
1425            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
1426            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
1427#if defined key_agrif
1428            IF ( .NOT.Agrif_Root() ) THEN
1429               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
1430            ENDIF
1431#endif
1432         ENDIF
1433         !
1434      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1435         !                                   ! -------------------
1436         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1437         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1438         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
1439         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
1440         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
1441         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
1442         !
1443         IF (.NOT.ln_bt_av) THEN
1444            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
1445            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
1446            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
1447            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
1448            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
1449            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
1450         ENDIF
1451#if defined key_agrif
1452         ! Save time integrated fluxes
1453         IF ( .NOT.Agrif_Root() ) THEN
1454            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
1455            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
1456         ENDIF
1457#endif
1458         IF( lwxios ) CALL iom_swap(      cxios_context          )
1459      ENDIF
1460      !
1461   END SUBROUTINE ts_rst
1462
1463
1464   SUBROUTINE dyn_spg_ts_init
1465      !!---------------------------------------------------------------------
1466      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1467      !!
1468      !! ** Purpose : Set time splitting options
1469      !!----------------------------------------------------------------------
1470      INTEGER  ::   ji ,jj              ! dummy loop indices
1471      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1472      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
1473      !!----------------------------------------------------------------------
1474      !
1475      ! Max courant number for ext. grav. waves
1476      !
1477      DO jj = 1, jpj
1478         DO ji =1, jpi
1479            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1480            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1481            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1482         END DO
1483      END DO
1484      !
1485      zcmax = MAXVAL( zcu(:,:) )
1486      CALL mpp_max( 'dynspg_ts', zcmax )
1487
1488      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1489      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1490     
1491      rdtbt = rdt / REAL( nn_baro , wp )
1492      zcmax = zcmax * rdtbt
1493      ! Print results
1494      IF(lwp) WRITE(numout,*)
1495      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1496      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1497      IF( ln_bt_auto ) THEN
1498         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
1499         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1500      ELSE
1501         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
1502      ENDIF
1503
1504      IF(ln_bt_av) THEN
1505         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
1506      ELSE
1507         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1508      ENDIF
1509      !
1510      !
1511      IF(ln_bt_fw) THEN
1512         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1513      ELSE
1514         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1515      ENDIF
1516      !
1517#if defined key_agrif
1518      ! Restrict the use of Agrif to the forward case only
1519!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1520#endif
1521      !
1522      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1523      SELECT CASE ( nn_bt_flt )
1524         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1525         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1526         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1527         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1528      END SELECT
1529      !
1530      IF(lwp) WRITE(numout,*) ' '
1531      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1532      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1533      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1534      !
1535      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1536      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1537         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1538      ENDIF
1539      !
1540      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1541         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1542      ENDIF
1543      IF( zcmax>0.9_wp ) THEN
1544         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1545      ENDIF
1546      !
1547      !                             ! Allocate time-splitting arrays
1548      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1549      !
1550      !                             ! read restart when needed
1551      CALL ts_rst( nit000, 'READ' )
1552      !
1553      IF( lwxios ) THEN
1554! define variables in restart file when writing with XIOS
1555         CALL iom_set_rstw_var_active('ub2_b')
1556         CALL iom_set_rstw_var_active('vb2_b')
1557         CALL iom_set_rstw_var_active('un_bf')
1558         CALL iom_set_rstw_var_active('vn_bf')
1559         !
1560         IF (.NOT.ln_bt_av) THEN
1561            CALL iom_set_rstw_var_active('sshbb_e')
1562            CALL iom_set_rstw_var_active('ubb_e')
1563            CALL iom_set_rstw_var_active('vbb_e')
1564            CALL iom_set_rstw_var_active('sshb_e')
1565            CALL iom_set_rstw_var_active('ub_e')
1566            CALL iom_set_rstw_var_active('vb_e')
1567         ENDIF
1568#if defined key_agrif
1569         ! Save time integrated fluxes
1570         IF ( .NOT.Agrif_Root() ) THEN
1571            CALL iom_set_rstw_var_active('ub2_i_b')
1572            CALL iom_set_rstw_var_active('vb2_i_b')
1573         ENDIF
1574#endif
1575      ENDIF
1576      !
1577   END SUBROUTINE dyn_spg_ts_init
1578
1579   !!======================================================================
1580END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.