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/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DYN/dynspg_ts.F90 @ 10852

Last change on this file since 10852 was 10852, checked in by smueller, 5 years ago

Renaming of subroutine sbc_tide of module sbctide to tide_update, transfer of this subroutine to module tide_mod, removal of the emptied module sbctide, and related adjustments (ticket #2194)

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