source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynspg_ts.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh

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