source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/OCE/DYN/dynspg_ts.F90 @ 11219

Last change on this file since 11219 was 11219, checked in by jchanut, 15 months ago

#2199
1) Define aditionnal arrays to correct the time interpolation of barotropic arrays in corners. Since multiple stages in the time interpolation are necessary, overlapping segments in corners give wrong results otherwise (corrects stage 2 in previous commit)..
2) Added subroutine to correct time extrapolated fluxes at bdy in time splitting routine (updates stage 3 in previous commit).
3) Completly remove non-specified open boundary case. Boundares are now exactly set from parent (no more filtering nor extrapolation for outgoing flows).
At this stage, use of nbondi, nbondj variables has been completly removed.

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