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

Last change on this file since 9528 was 9528, checked in by gm, 2 years ago

dev_merge_2017: add to new vorticity scheme (EET= EEN using e3t instead of e3f, and ENT= energy conserving at T-point). NB: not used in ref config.

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