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

source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 10481

Last change on this file since 10481 was 10481, checked in by mathiot, 5 years ago

fix ticket #2200

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