source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90 @ 11541

Last change on this file since 11541 was 11541, checked in by mathiot, 12 months ago

ENHANCE-02_ISF: simplify use of ln_isf, add extra comments + minor changes (ticket #2142)

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