source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9554

Last change on this file since 9554 was 9554, checked in by mathiot, 2 years ago

Correction bug with isf and flux form (mask issue) + change mask used for ubaro output

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