source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90 @ 11381

Last change on this file since 11381 was 11381, checked in by girrmann, 2 years ago

dev_r10984_HPC-13 : repair faulty commit [11380], see #2308

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