source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/dynspg_ts.F90 @ 11053

Last change on this file since 11053 was 11053, checked in by davestorkey, 19 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Merge in latest changes from main branch and finish conversion of "h" variables. NB. This version still doesn't work!

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