source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN/dynspg_ts.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

  • Property svn:keywords set to Id
File size: 77.0 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      IF( lk_mpp                )   CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
115      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('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         !                                                !* Update the forcing (BDY and tides)
714         !                                                !  ------------------
715         ! Update only tidal forcing at open boundaries
716         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
717         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
718         !
719         ! Set extrapolation coefficients for predictor step:
720         IF ((jn<3).AND.ll_init) THEN      ! Forward           
721           za1 = 1._wp                                         
722           za2 = 0._wp                       
723           za3 = 0._wp                       
724         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
725           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
726           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
727           za3 =  0.281105_wp              ! za3 = bet
728         ENDIF
729
730         ! Extrapolate barotropic velocities at step jit+0.5:
731         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
732         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
733
734         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
735            !                                             !  ------------------
736            ! Extrapolate Sea Level at step jit+0.5:
737            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
738           
739            ! set wetting & drying mask at tracer points for this barotropic sub-step
740            IF ( ln_wd_dl ) THEN 
741               !
742               IF ( ln_wd_dl_rmp ) THEN
743                  DO jj = 1, jpj                                 
744                     DO ji = 1, jpi   ! vector opt. 
745                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
746!                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN
747                           ztwdmask(ji,jj) = 1._wp
748                        ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
749                           ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) ) 
750                        ELSE
751                           ztwdmask(ji,jj) = 0._wp
752                        END IF
753                     END DO
754                  END DO
755               ELSE
756                  DO jj = 1, jpj                                 
757                     DO ji = 1, jpi   ! vector opt. 
758                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
759                           ztwdmask(ji,jj) = 1._wp
760                        ELSE
761                           ztwdmask(ji,jj) = 0._wp
762                        ENDIF
763                     END DO
764                  END DO
765               ENDIF 
766               !
767            ENDIF 
768            !
769            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
770               DO ji = 2, fs_jpim1   ! Vector opt.
771                  zwx(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
772                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
773                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
774                  zwy(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
775                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
776                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
777               END DO
778            END DO
779            CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp )
780            !
781            zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
782            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:)
783            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
784         ELSE
785            zhup2_e(:,:) = hu_n(:,:)
786            zhvp2_e(:,:) = hv_n(:,:)
787            zhtp2_e(:,:) = ht_n(:,:)
788         ENDIF
789         !                                                !* after ssh
790         !                                                !  -----------
791         ! One should enforce volume conservation at open boundaries here
792         ! considering fluxes below:
793         !
794         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
795         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
796         !
797#if defined key_agrif
798         ! Set fluxes during predictor step to ensure volume conservation
799         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
800            IF((nbondi == -1).OR.(nbondi == 2)) THEN
801               DO jj = 1, jpj
802                  zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
803                  zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
804               END DO
805            ENDIF
806            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
807               DO jj=1,jpj
808                  zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
809                  zwy(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj)
810               END DO
811            ENDIF
812            IF((nbondj == -1).OR.(nbondj == 2)) THEN
813               DO ji=1,jpi
814                  zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
815                  zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
816               END DO
817            ENDIF
818            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
819               DO ji=1,jpi
820                  zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
821                  zwx(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1)
822               END DO
823            ENDIF
824         ENDIF
825#endif
826         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
827
828         IF ( ln_wd_dl ) THEN 
829            !
830            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
831            !
832            DO jj = 1, jpjm1                                 
833               DO ji = 1, jpim1   
834                  IF ( zwx(ji,jj) > 0.0 ) THEN
835                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj) 
836                  ELSE
837                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 
838                  END IF
839                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
840                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
841
842                  IF ( zwy(ji,jj) > 0.0 ) THEN
843                     zvwdmask(ji, jj) = ztwdmask(ji, jj  )
844                  ELSE
845                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 
846                  END IF
847                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 
848                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
849               END DO
850            END DO
851            !
852         ENDIF   
853         
854         ! Sum over sub-time-steps to compute advective velocities
855         za2 = wgtbtp2(jn)
856         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
857         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
858         
859         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
860         IF ( ln_wd_dl_bc ) THEN
861            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
862            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
863         END IF 
864
865         ! Set next sea level:
866         DO jj = 2, jpjm1                                 
867            DO ji = fs_2, fs_jpim1   ! vector opt.
868               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
869                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
870            END DO
871         END DO
872         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
873         
874         CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T',  1._wp )
875
876         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
877         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
878#if defined key_agrif
879         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
880#endif
881         
882         ! Sea Surface Height at u-,v-points (vvl case only)
883         IF( .NOT.ln_linssh ) THEN                               
884            DO jj = 2, jpjm1
885               DO ji = 2, jpim1      ! NO Vector Opt.
886                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
887                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
888                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
889                  zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
890                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
891                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
892               END DO
893            END DO
894            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
895         ENDIF   
896         !                                 
897         ! Half-step back interpolation of SSH for surface pressure computation:
898         !----------------------------------------------------------------------
899         IF ((jn==1).AND.ll_init) THEN
900           za0=1._wp                        ! Forward-backward
901           za1=0._wp                           
902           za2=0._wp
903           za3=0._wp
904         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
905           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
906           za1=-0.1666666666666_wp          ! za1 = gam
907           za2= 0.0833333333333_wp          ! za2 = eps
908           za3= 0._wp             
909         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
910            IF (rn_bt_alpha==0._wp) THEN
911               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps
912               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps
913               za2=0.088_wp                 ! za2 = gam
914               za3=0.013_wp                 ! za3 = eps
915            ELSE
916               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
917               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
918               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
919               za1 = 1._wp - za0 - zgamma - zepsilon
920               za2 = zgamma
921               za3 = zepsilon
922            ENDIF
923         ENDIF
924         !
925         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
926            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
927         
928         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters
929           DO jj = 2, jpjm1
930              DO ji = 2, jpim1 
931                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
932                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            &
933                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
934                     &                                                             > rn_wdmin1 + rn_wdmin2
935                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
936                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
937                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
938   
939                IF(ll_tmp1) THEN
940                  zcpx(ji,jj) = 1.0_wp
941                ELSE IF(ll_tmp2) THEN
942                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here
943                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
944                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) )
945                ELSE
946                  zcpx(ji,jj) = 0._wp
947                ENDIF
948                !
949                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
950                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            &
951                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
952                     &                                                             > rn_wdmin1 + rn_wdmin2
953                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
954                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
955                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
956   
957                IF(ll_tmp1) THEN
958                  zcpy(ji,jj) = 1.0_wp
959                ELSEIF(ll_tmp2) THEN
960                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
961                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
962                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) )
963                ELSE
964                  zcpy(ji,jj) = 0._wp
965                ENDIF
966              END DO
967           END DO
968         ENDIF
969         !
970         ! Compute associated depths at U and V points:
971         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form
972            !                                       
973            DO jj = 2, jpjm1                           
974               DO ji = 2, jpim1
975                  zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
976                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
977                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
978                  zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
979                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
980                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
981                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
982                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
983               END DO
984            END DO
985            !
986         ENDIF
987         !
988         ! Add Coriolis trend:
989         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
990         ! at each time step. We however keep them constant here for optimization.
991         ! Recall that zwx and zwy arrays hold fluxes at this stage:
992         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
993         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
994         !
995         SELECT CASE( nvor_scheme )
996         CASE( np_ENT )             ! energy conserving scheme (t-point)
997         DO jj = 2, jpjm1
998            DO ji = 2, jpim1   ! vector opt.
999
1000               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
1001               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
1002           
1003               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   &
1004                  &               * (  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) )   &
1005                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   )
1006                  !
1007               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1008                  &               * (  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) )   & 
1009                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   ) 
1010            END DO 
1011         END DO 
1012         !         
1013         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point)
1014            DO jj = 2, jpjm1
1015               DO ji = fs_2, fs_jpim1   ! vector opt.
1016                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
1017                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
1018                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
1019                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
1020                  zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1021                  zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1022               END DO
1023            END DO
1024            !
1025         CASE( np_ENS )             ! enstrophy conserving scheme (f-point)
1026            DO jj = 2, jpjm1
1027               DO ji = fs_2, fs_jpim1   ! vector opt.
1028                  zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
1029                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
1030                  zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
1031                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
1032                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1033                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1034               END DO
1035            END DO
1036            !
1037         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f)
1038            DO jj = 2, jpjm1
1039               DO ji = fs_2, fs_jpim1   ! vector opt.
1040                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  &
1041                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  &
1042                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  & 
1043                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  )
1044                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  & 
1045                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  &
1046                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  & 
1047                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  )
1048               END DO
1049            END DO
1050            !
1051         END SELECT
1052         !
1053         ! Add tidal astronomical forcing if defined
1054         IF ( ln_tide .AND. ln_tide_pot ) THEN
1055            DO jj = 2, jpjm1
1056               DO ji = fs_2, fs_jpim1   ! vector opt.
1057                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1058                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1059                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1060                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1061               END DO
1062            END DO
1063         ENDIF
1064         !
1065         ! Add bottom stresses:
1066!jth do implicitly instead
1067         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1068            DO jj = 2, jpjm1
1069               DO ji = fs_2, fs_jpim1   ! vector opt.
1070                  zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1071                  zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1072               END DO
1073            END DO
1074         ENDIF 
1075         !
1076         ! Surface pressure trend:
1077         IF( ln_wd_il ) THEN
1078           DO jj = 2, jpjm1
1079              DO ji = 2, jpim1 
1080                 ! Add surface pressure gradient
1081                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1082                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1083                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 
1084                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1085              END DO
1086           END DO
1087         ELSE
1088           DO jj = 2, jpjm1
1089              DO ji = fs_2, fs_jpim1   ! vector opt.
1090                 ! Add surface pressure gradient
1091                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1092                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1093                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1094                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1095              END DO
1096           END DO
1097         END IF
1098
1099         !
1100         ! Set next velocities:
1101         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
1102            DO jj = 2, jpjm1
1103               DO ji = fs_2, fs_jpim1   ! vector opt.
1104                  ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
1105                            &     + rdtbt * (                      zwx(ji,jj)   &
1106                            &                                 + zu_trd(ji,jj)   &
1107                            &                                 + zu_frc(ji,jj) ) & 
1108                            &   ) * ssumask(ji,jj)
1109
1110                  va_e(ji,jj) = (                                 vn_e(ji,jj)   &
1111                            &     + rdtbt * (                      zwy(ji,jj)   &
1112                            &                                 + zv_trd(ji,jj)   &
1113                            &                                 + zv_frc(ji,jj) ) &
1114                            &   ) * ssvmask(ji,jj)
1115 
1116               END DO
1117            END DO
1118            !
1119         ELSE                           !* Flux form
1120            DO jj = 2, jpjm1
1121               DO ji = fs_2, fs_jpim1   ! vector opt.
1122
1123                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1124                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1125
1126                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1127                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1128
1129                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
1130                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
1131                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
1132                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
1133                            &   ) * zhura
1134
1135                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
1136                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
1137                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
1138                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
1139                            &   ) * zhvra
1140               END DO
1141            END DO
1142         ENDIF
1143!jth implicit bottom friction:
1144         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1145            DO jj = 2, jpjm1
1146               DO ji = fs_2, fs_jpim1   ! vector opt.
1147                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1148                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1149               END DO
1150            END DO
1151         ENDIF
1152
1153         
1154         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
1155            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1156            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1157            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1158            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1159            !
1160         ENDIF
1161         !                                             !* domain lateral boundary
1162         CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1163         !
1164         !                                                 ! open boundaries
1165         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1166#if defined key_agrif                                                           
1167         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
1168#endif
1169         !                                             !* Swap
1170         !                                             !  ----
1171         ubb_e  (:,:) = ub_e  (:,:)
1172         ub_e   (:,:) = un_e  (:,:)
1173         un_e   (:,:) = ua_e  (:,:)
1174         !
1175         vbb_e  (:,:) = vb_e  (:,:)
1176         vb_e   (:,:) = vn_e  (:,:)
1177         vn_e   (:,:) = va_e  (:,:)
1178         !
1179         sshbb_e(:,:) = sshb_e(:,:)
1180         sshb_e (:,:) = sshn_e(:,:)
1181         sshn_e (:,:) = ssha_e(:,:)
1182
1183         !                                             !* Sum over whole bt loop
1184         !                                             !  ----------------------
1185         za1 = wgtbtp1(jn)                                   
1186         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
1187            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
1188            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
1189         ELSE                                       ! Sum transports
1190            IF ( .NOT.ln_wd_dl ) THEN 
1191               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
1192               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
1193            ELSE
1194               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
1195               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
1196            END IF
1197         ENDIF
1198         !                                          ! Sum sea level
1199         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1200
1201         !                                                 ! ==================== !
1202      END DO                                               !        end loop      !
1203      !                                                    ! ==================== !
1204      ! -----------------------------------------------------------------------------
1205      ! Phase 3. update the general trend with the barotropic trend
1206      ! -----------------------------------------------------------------------------
1207      !
1208      ! Set advection velocity correction:
1209      IF (ln_bt_fw) THEN
1210         zwx(:,:) = un_adv(:,:)
1211         zwy(:,:) = vn_adv(:,:)
1212         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1213            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1214            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1215            !
1216            ! Update corrective fluxes for next time step:
1217            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1218            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1219         ELSE
1220            un_bf(:,:) = 0._wp
1221            vn_bf(:,:) = 0._wp 
1222         END IF         
1223         ! Save integrated transport for next computation
1224         ub2_b(:,:) = zwx(:,:)
1225         vb2_b(:,:) = zwy(:,:)
1226      ENDIF
1227
1228
1229      !
1230      ! Update barotropic trend:
1231      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1232         DO jk=1,jpkm1
1233            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1234            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1235         END DO
1236      ELSE
1237         ! At this stage, ssha has been corrected: compute new depths at velocity points
1238         DO jj = 1, jpjm1
1239            DO ji = 1, jpim1      ! NO Vector Opt.
1240               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
1241                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      &
1242                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1243               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
1244                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
1245                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1246            END DO
1247         END DO
1248         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1249         !
1250         DO jk=1,jpkm1
1251            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1252            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1253         END DO
1254         ! Save barotropic velocities not transport:
1255         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1256         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1257      ENDIF
1258
1259
1260      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
1261      DO jk = 1, jpkm1
1262         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)
1263         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1264      END DO
1265
1266      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1267         DO jk = 1, jpkm1
1268            un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1269                       & + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk) 
1270            vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) & 
1271                       & + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk) 
1272         END DO
1273      END IF
1274
1275     
1276      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current
1277      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current
1278      !
1279#if defined key_agrif
1280      ! Save time integrated fluxes during child grid integration
1281      ! (used to update coarse grid transports at next time step)
1282      !
1283      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1284         IF( Agrif_NbStepint() == 0 ) THEN
1285            ub2_i_b(:,:) = 0._wp
1286            vb2_i_b(:,:) = 0._wp
1287         END IF
1288         !
1289         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1290         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1291         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1292      ENDIF
1293#endif     
1294      !                                   !* write time-spliting arrays in the restart
1295      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
1296      !
1297      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
1298      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1299      !
1300      IF( ln_diatmb ) THEN
1301         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
1302         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
1303      ENDIF
1304      !
1305   END SUBROUTINE dyn_spg_ts
1306
1307
1308   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1309      !!---------------------------------------------------------------------
1310      !!                   ***  ROUTINE ts_wgt  ***
1311      !!
1312      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1313      !!----------------------------------------------------------------------
1314      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
1315      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
1316      INTEGER, INTENT(inout) :: jpit      ! cycle length   
1317      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
1318                                                         zwgt2    ! Secondary weights
1319     
1320      INTEGER ::  jic, jn, ji                      ! temporary integers
1321      REAL(wp) :: za1, za2
1322      !!----------------------------------------------------------------------
1323
1324      zwgt1(:) = 0._wp
1325      zwgt2(:) = 0._wp
1326
1327      ! Set time index when averaged value is requested
1328      IF (ll_fw) THEN
1329         jic = nn_baro
1330      ELSE
1331         jic = 2 * nn_baro
1332      ENDIF
1333
1334      ! Set primary weights:
1335      IF (ll_av) THEN
1336           ! Define simple boxcar window for primary weights
1337           ! (width = nn_baro, centered around jic)     
1338         SELECT CASE ( nn_bt_flt )
1339              CASE( 0 )  ! No averaging
1340                 zwgt1(jic) = 1._wp
1341                 jpit = jic
1342
1343              CASE( 1 )  ! Boxcar, width = nn_baro
1344                 DO jn = 1, 3*nn_baro
1345                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1346                    IF (za1 < 0.5_wp) THEN
1347                      zwgt1(jn) = 1._wp
1348                      jpit = jn
1349                    ENDIF
1350                 ENDDO
1351
1352              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
1353                 DO jn = 1, 3*nn_baro
1354                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1355                    IF (za1 < 1._wp) THEN
1356                      zwgt1(jn) = 1._wp
1357                      jpit = jn
1358                    ENDIF
1359                 ENDDO
1360              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1361         END SELECT
1362
1363      ELSE ! No time averaging
1364         zwgt1(jic) = 1._wp
1365         jpit = jic
1366      ENDIF
1367   
1368      ! Set secondary weights
1369      DO jn = 1, jpit
1370        DO ji = jn, jpit
1371             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1372        END DO
1373      END DO
1374
1375      ! Normalize weigths:
1376      za1 = 1._wp / SUM(zwgt1(1:jpit))
1377      za2 = 1._wp / SUM(zwgt2(1:jpit))
1378      DO jn = 1, jpit
1379        zwgt1(jn) = zwgt1(jn) * za1
1380        zwgt2(jn) = zwgt2(jn) * za2
1381      END DO
1382      !
1383   END SUBROUTINE ts_wgt
1384
1385
1386   SUBROUTINE ts_rst( kt, cdrw )
1387      !!---------------------------------------------------------------------
1388      !!                   ***  ROUTINE ts_rst  ***
1389      !!
1390      !! ** Purpose : Read or write time-splitting arrays in restart file
1391      !!----------------------------------------------------------------------
1392      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1393      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1394      !!----------------------------------------------------------------------
1395      !
1396      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1397         !                                   ! ---------------
1398         IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN    !* Read the restart file
1399            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
1400            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
1401            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
1402            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
1403            IF( .NOT.ln_bt_av ) THEN
1404               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
1405               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
1406               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
1407               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
1408               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
1409               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
1410            ENDIF
1411#if defined key_agrif
1412            ! Read time integrated fluxes
1413            IF ( .NOT.Agrif_Root() ) THEN
1414               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
1415               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
1416            ENDIF
1417#endif
1418         ELSE                                   !* Start from rest
1419            IF(lwp) WRITE(numout,*)
1420            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
1421            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
1422            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
1423            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
1424#if defined key_agrif
1425            IF ( .NOT.Agrif_Root() ) THEN
1426               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
1427            ENDIF
1428#endif
1429         ENDIF
1430         !
1431      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1432         !                                   ! -------------------
1433         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1434         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1435         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
1436         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
1437         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
1438         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
1439         !
1440         IF (.NOT.ln_bt_av) THEN
1441            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
1442            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
1443            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
1444            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
1445            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
1446            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
1447         ENDIF
1448#if defined key_agrif
1449         ! Save time integrated fluxes
1450         IF ( .NOT.Agrif_Root() ) THEN
1451            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
1452            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
1453         ENDIF
1454#endif
1455         IF( lwxios ) CALL iom_swap(      cxios_context          )
1456      ENDIF
1457      !
1458   END SUBROUTINE ts_rst
1459
1460
1461   SUBROUTINE dyn_spg_ts_init
1462      !!---------------------------------------------------------------------
1463      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1464      !!
1465      !! ** Purpose : Set time splitting options
1466      !!----------------------------------------------------------------------
1467      INTEGER  ::   ji ,jj              ! dummy loop indices
1468      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1469      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
1470      !!----------------------------------------------------------------------
1471      !
1472      ! Max courant number for ext. grav. waves
1473      !
1474      DO jj = 1, jpj
1475         DO ji =1, jpi
1476            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1477            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1478            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1479         END DO
1480      END DO
1481      !
1482      zcmax = MAXVAL( zcu(:,:) )
1483      IF( lk_mpp )   CALL mpp_max( 'dynspg_ts', zcmax )
1484
1485      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1486      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1487     
1488      rdtbt = rdt / REAL( nn_baro , wp )
1489      zcmax = zcmax * rdtbt
1490      ! Print results
1491      IF(lwp) WRITE(numout,*)
1492      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1493      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1494      IF( ln_bt_auto ) THEN
1495         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
1496         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1497      ELSE
1498         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
1499      ENDIF
1500
1501      IF(ln_bt_av) THEN
1502         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
1503      ELSE
1504         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1505      ENDIF
1506      !
1507      !
1508      IF(ln_bt_fw) THEN
1509         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1510      ELSE
1511         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1512      ENDIF
1513      !
1514#if defined key_agrif
1515      ! Restrict the use of Agrif to the forward case only
1516!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1517#endif
1518      !
1519      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1520      SELECT CASE ( nn_bt_flt )
1521         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1522         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1523         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1524         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1525      END SELECT
1526      !
1527      IF(lwp) WRITE(numout,*) ' '
1528      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1529      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1530      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1531      !
1532      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1533      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1534         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1535      ENDIF
1536      !
1537      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1538         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1539      ENDIF
1540      IF( zcmax>0.9_wp ) THEN
1541         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1542      ENDIF
1543      !
1544      !                             ! Allocate time-splitting arrays
1545      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1546      !
1547      !                             ! read restart when needed
1548      CALL ts_rst( nit000, 'READ' )
1549      !
1550      IF( lwxios ) THEN
1551! define variables in restart file when writing with XIOS
1552         CALL iom_set_rstw_var_active('ub2_b')
1553         CALL iom_set_rstw_var_active('vb2_b')
1554         CALL iom_set_rstw_var_active('un_bf')
1555         CALL iom_set_rstw_var_active('vn_bf')
1556         !
1557         IF (.NOT.ln_bt_av) THEN
1558            CALL iom_set_rstw_var_active('sshbb_e')
1559            CALL iom_set_rstw_var_active('ubb_e')
1560            CALL iom_set_rstw_var_active('vbb_e')
1561            CALL iom_set_rstw_var_active('sshb_e')
1562            CALL iom_set_rstw_var_active('ub_e')
1563            CALL iom_set_rstw_var_active('vb_e')
1564         ENDIF
1565#if defined key_agrif
1566         ! Save time integrated fluxes
1567         IF ( .NOT.Agrif_Root() ) THEN
1568            CALL iom_set_rstw_var_active('ub2_i_b')
1569            CALL iom_set_rstw_var_active('vb2_i_b')
1570         ENDIF
1571#endif
1572      ENDIF
1573      !
1574   END SUBROUTINE dyn_spg_ts_init
1575
1576   !!======================================================================
1577END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.