source: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg_ts.F90 @ 11613

Last change on this file since 11613 was 11613, checked in by davestorkey, 14 months ago

UKMO/NEMO_4.0_momentum_trends branch : first set of code changes.

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