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/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN – NEMO

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

Last change on this file since 9939 was 9939, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

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