New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in NEMO/branches/UKMO/dev_r9885_proto_GO8_package/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/dev_r9885_proto_GO8_package/src/OCE/DYN/dynspg_ts.F90 @ 9892

Last change on this file since 9892 was 9892, checked in by davestorkey, 6 years ago

UKMO/dev_r9885_proto_GO8_package branch: clear SVN keywords

File size: 76.7 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 licence     (./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( 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( 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( 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( 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( 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( 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!jth implicit bottom friction:
1117                  IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1118                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1119                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1120                  ENDIF
1121
1122               END DO
1123            END DO
1124            !
1125         ELSE                           !* Flux form
1126            DO jj = 2, jpjm1
1127               DO ji = fs_2, fs_jpim1   ! vector opt.
1128
1129                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1130                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1131
1132                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1133                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1134
1135                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
1136                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
1137                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
1138                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
1139                            &   ) * zhura
1140
1141                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
1142                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
1143                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
1144                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
1145                            &   ) * zhvra
1146               END DO
1147            END DO
1148         ENDIF
1149
1150         
1151         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
1152            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1153            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1154            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1155            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1156            !
1157         ENDIF
1158         !                                             !* domain lateral boundary
1159         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1160         !
1161         !                                                 ! open boundaries
1162         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1163#if defined key_agrif                                                           
1164         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
1165#endif
1166         !                                             !* Swap
1167         !                                             !  ----
1168         ubb_e  (:,:) = ub_e  (:,:)
1169         ub_e   (:,:) = un_e  (:,:)
1170         un_e   (:,:) = ua_e  (:,:)
1171         !
1172         vbb_e  (:,:) = vb_e  (:,:)
1173         vb_e   (:,:) = vn_e  (:,:)
1174         vn_e   (:,:) = va_e  (:,:)
1175         !
1176         sshbb_e(:,:) = sshb_e(:,:)
1177         sshb_e (:,:) = sshn_e(:,:)
1178         sshn_e (:,:) = ssha_e(:,:)
1179
1180         !                                             !* Sum over whole bt loop
1181         !                                             !  ----------------------
1182         za1 = wgtbtp1(jn)                                   
1183         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
1184            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
1185            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
1186         ELSE                                       ! Sum transports
1187            IF ( .NOT.ln_wd_dl ) THEN 
1188               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
1189               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
1190            ELSE
1191               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
1192               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
1193            END IF
1194         ENDIF
1195         !                                          ! Sum sea level
1196         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1197
1198         !                                                 ! ==================== !
1199      END DO                                               !        end loop      !
1200      !                                                    ! ==================== !
1201      ! -----------------------------------------------------------------------------
1202      ! Phase 3. update the general trend with the barotropic trend
1203      ! -----------------------------------------------------------------------------
1204      !
1205      ! Set advection velocity correction:
1206      IF (ln_bt_fw) THEN
1207         zwx(:,:) = un_adv(:,:)
1208         zwy(:,:) = vn_adv(:,:)
1209         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1210            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1211            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1212            !
1213            ! Update corrective fluxes for next time step:
1214            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1215            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1216         ELSE
1217            un_bf(:,:) = 0._wp
1218            vn_bf(:,:) = 0._wp 
1219         END IF         
1220         ! Save integrated transport for next computation
1221         ub2_b(:,:) = zwx(:,:)
1222         vb2_b(:,:) = zwy(:,:)
1223      ENDIF
1224
1225
1226      !
1227      ! Update barotropic trend:
1228      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1229         DO jk=1,jpkm1
1230            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1231            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1232         END DO
1233      ELSE
1234         ! At this stage, ssha has been corrected: compute new depths at velocity points
1235         DO jj = 1, jpjm1
1236            DO ji = 1, jpim1      ! NO Vector Opt.
1237               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
1238                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      &
1239                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1240               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
1241                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
1242                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1243            END DO
1244         END DO
1245         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1246         !
1247         DO jk=1,jpkm1
1248            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1249            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1250         END DO
1251         ! Save barotropic velocities not transport:
1252         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1253         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1254      ENDIF
1255
1256
1257      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
1258      DO jk = 1, jpkm1
1259         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)
1260         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1261      END DO
1262
1263      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1264         DO jk = 1, jpkm1
1265            un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1266                       & + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk) 
1267            vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) & 
1268                       & + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk) 
1269         END DO
1270      END IF
1271
1272     
1273      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current
1274      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current
1275      !
1276#if defined key_agrif
1277      ! Save time integrated fluxes during child grid integration
1278      ! (used to update coarse grid transports at next time step)
1279      !
1280      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1281         IF( Agrif_NbStepint() == 0 ) THEN
1282            ub2_i_b(:,:) = 0._wp
1283            vb2_i_b(:,:) = 0._wp
1284         END IF
1285         !
1286         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1287         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1288         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1289      ENDIF
1290#endif     
1291      !                                   !* write time-spliting arrays in the restart
1292      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
1293      !
1294      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
1295      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1296      !
1297      IF( ln_diatmb ) THEN
1298         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
1299         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
1300      ENDIF
1301      !
1302   END SUBROUTINE dyn_spg_ts
1303
1304
1305   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1306      !!---------------------------------------------------------------------
1307      !!                   ***  ROUTINE ts_wgt  ***
1308      !!
1309      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1310      !!----------------------------------------------------------------------
1311      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
1312      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
1313      INTEGER, INTENT(inout) :: jpit      ! cycle length   
1314      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
1315                                                         zwgt2    ! Secondary weights
1316     
1317      INTEGER ::  jic, jn, ji                      ! temporary integers
1318      REAL(wp) :: za1, za2
1319      !!----------------------------------------------------------------------
1320
1321      zwgt1(:) = 0._wp
1322      zwgt2(:) = 0._wp
1323
1324      ! Set time index when averaged value is requested
1325      IF (ll_fw) THEN
1326         jic = nn_baro
1327      ELSE
1328         jic = 2 * nn_baro
1329      ENDIF
1330
1331      ! Set primary weights:
1332      IF (ll_av) THEN
1333           ! Define simple boxcar window for primary weights
1334           ! (width = nn_baro, centered around jic)     
1335         SELECT CASE ( nn_bt_flt )
1336              CASE( 0 )  ! No averaging
1337                 zwgt1(jic) = 1._wp
1338                 jpit = jic
1339
1340              CASE( 1 )  ! Boxcar, width = nn_baro
1341                 DO jn = 1, 3*nn_baro
1342                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1343                    IF (za1 < 0.5_wp) THEN
1344                      zwgt1(jn) = 1._wp
1345                      jpit = jn
1346                    ENDIF
1347                 ENDDO
1348
1349              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
1350                 DO jn = 1, 3*nn_baro
1351                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1352                    IF (za1 < 1._wp) THEN
1353                      zwgt1(jn) = 1._wp
1354                      jpit = jn
1355                    ENDIF
1356                 ENDDO
1357              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1358         END SELECT
1359
1360      ELSE ! No time averaging
1361         zwgt1(jic) = 1._wp
1362         jpit = jic
1363      ENDIF
1364   
1365      ! Set secondary weights
1366      DO jn = 1, jpit
1367        DO ji = jn, jpit
1368             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1369        END DO
1370      END DO
1371
1372      ! Normalize weigths:
1373      za1 = 1._wp / SUM(zwgt1(1:jpit))
1374      za2 = 1._wp / SUM(zwgt2(1:jpit))
1375      DO jn = 1, jpit
1376        zwgt1(jn) = zwgt1(jn) * za1
1377        zwgt2(jn) = zwgt2(jn) * za2
1378      END DO
1379      !
1380   END SUBROUTINE ts_wgt
1381
1382
1383   SUBROUTINE ts_rst( kt, cdrw )
1384      !!---------------------------------------------------------------------
1385      !!                   ***  ROUTINE ts_rst  ***
1386      !!
1387      !! ** Purpose : Read or write time-splitting arrays in restart file
1388      !!----------------------------------------------------------------------
1389      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1390      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1391      !!----------------------------------------------------------------------
1392      !
1393      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1394         !                                   ! ---------------
1395         IF( ln_rstart .AND. ln_bt_fw ) THEN    !* Read the restart file
1396            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
1397            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
1398            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
1399            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
1400            IF( .NOT.ln_bt_av ) THEN
1401               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
1402               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
1403               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
1404               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
1405               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
1406               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
1407            ENDIF
1408#if defined key_agrif
1409            ! Read time integrated fluxes
1410            IF ( .NOT.Agrif_Root() ) THEN
1411               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
1412               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
1413            ENDIF
1414#endif
1415         ELSE                                   !* Start from rest
1416            IF(lwp) WRITE(numout,*)
1417            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
1418            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
1419            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
1420            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
1421#if defined key_agrif
1422            IF ( .NOT.Agrif_Root() ) THEN
1423               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
1424            ENDIF
1425#endif
1426         ENDIF
1427         !
1428      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1429         !                                   ! -------------------
1430         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1431         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1432         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
1433         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
1434         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
1435         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
1436         !
1437         IF (.NOT.ln_bt_av) THEN
1438            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
1439            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
1440            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
1441            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
1442            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
1443            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
1444         ENDIF
1445#if defined key_agrif
1446         ! Save time integrated fluxes
1447         IF ( .NOT.Agrif_Root() ) THEN
1448            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
1449            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
1450         ENDIF
1451#endif
1452         IF( lwxios ) CALL iom_swap(      cxios_context          )
1453      ENDIF
1454      !
1455   END SUBROUTINE ts_rst
1456
1457
1458   SUBROUTINE dyn_spg_ts_init
1459      !!---------------------------------------------------------------------
1460      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1461      !!
1462      !! ** Purpose : Set time splitting options
1463      !!----------------------------------------------------------------------
1464      INTEGER  ::   ji ,jj              ! dummy loop indices
1465      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1466      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
1467      !!----------------------------------------------------------------------
1468      !
1469      ! Max courant number for ext. grav. waves
1470      !
1471      DO jj = 1, jpj
1472         DO ji =1, jpi
1473            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1474            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1475            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1476         END DO
1477      END DO
1478      !
1479      zcmax = MAXVAL( zcu(:,:) )
1480      IF( lk_mpp )   CALL mpp_max( zcmax )
1481
1482      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1483      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1484     
1485      rdtbt = rdt / REAL( nn_baro , wp )
1486      zcmax = zcmax * rdtbt
1487      ! Print results
1488      IF(lwp) WRITE(numout,*)
1489      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1490      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1491      IF( ln_bt_auto ) THEN
1492         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
1493         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1494      ELSE
1495         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
1496      ENDIF
1497
1498      IF(ln_bt_av) THEN
1499         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
1500      ELSE
1501         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1502      ENDIF
1503      !
1504      !
1505      IF(ln_bt_fw) THEN
1506         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1507      ELSE
1508         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1509      ENDIF
1510      !
1511#if defined key_agrif
1512      ! Restrict the use of Agrif to the forward case only
1513!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1514#endif
1515      !
1516      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1517      SELECT CASE ( nn_bt_flt )
1518         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1519         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1520         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1521         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1522      END SELECT
1523      !
1524      IF(lwp) WRITE(numout,*) ' '
1525      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1526      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1527      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1528      !
1529      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1530      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1531         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1532      ENDIF
1533      !
1534      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1535         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1536      ENDIF
1537      IF( zcmax>0.9_wp ) THEN
1538         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1539      ENDIF
1540      !
1541      !                             ! Allocate time-splitting arrays
1542      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1543      !
1544      !                             ! read restart when needed
1545      CALL ts_rst( nit000, 'READ' )
1546      !
1547      IF( lwxios ) THEN
1548! define variables in restart file when writing with XIOS
1549         CALL iom_set_rstw_var_active('ub2_b')
1550         CALL iom_set_rstw_var_active('vb2_b')
1551         CALL iom_set_rstw_var_active('un_bf')
1552         CALL iom_set_rstw_var_active('vn_bf')
1553         !
1554         IF (.NOT.ln_bt_av) THEN
1555            CALL iom_set_rstw_var_active('sshbb_e')
1556            CALL iom_set_rstw_var_active('ubb_e')
1557            CALL iom_set_rstw_var_active('vbb_e')
1558            CALL iom_set_rstw_var_active('sshb_e')
1559            CALL iom_set_rstw_var_active('ub_e')
1560            CALL iom_set_rstw_var_active('vb_e')
1561         ENDIF
1562#if defined key_agrif
1563         ! Save time integrated fluxes
1564         IF ( .NOT.Agrif_Root() ) THEN
1565            CALL iom_set_rstw_var_active('ub2_i_b')
1566            CALL iom_set_rstw_var_active('vb2_i_b')
1567         ENDIF
1568#endif
1569      ENDIF
1570      !
1571   END SUBROUTINE dyn_spg_ts_init
1572
1573   !!======================================================================
1574END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.