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

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

#2199
1) Make sponge independent of sub-domain size. Partially masked open boundary segments are not taken into account anymore. To do so, sponge coefficients should be read in a file for realistic applications (then nesting tools need to be modified accordingly).
2) Replace East-West-North-South barotropic data arrays by a global 2d array. Then apply barotropic open boundary conditions thanks to mi0/mi1, mj0/mj1 indexes.
3) Call AGRIF bdy update one more time in dynspg_ts during extrapolation phase. This removes a dozen lines of code in dynspg_ts routine.

  • 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#if  defined key_agrif
800         ! Set fluxes during predictor step to ensure volume conservation
801         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts( jn )
802#endif     
803         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
804         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
805         !
806         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
807
808         IF ( ln_wd_dl ) THEN 
809            !
810            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
811            !
812            DO jj = 1, jpjm1                                 
813               DO ji = 1, jpim1   
814                  IF ( zwx(ji,jj) > 0.0 ) THEN
815                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj) 
816                  ELSE
817                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 
818                  END IF
819                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
820                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
821
822                  IF ( zwy(ji,jj) > 0.0 ) THEN
823                     zvwdmask(ji, jj) = ztwdmask(ji, jj  )
824                  ELSE
825                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 
826                  END IF
827                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 
828                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
829               END DO
830            END DO
831            !
832         ENDIF   
833         
834         ! Sum over sub-time-steps to compute advective velocities
835         za2 = wgtbtp2(jn)
836         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
837         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
838         
839         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
840         IF ( ln_wd_dl_bc ) THEN
841            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
842            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
843         END IF 
844
845         ! Set next sea level:
846         DO jj = 2, jpjm1                                 
847            DO ji = fs_2, fs_jpim1   ! vector opt.
848               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
849                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
850            END DO
851         END DO
852         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
853         
854         CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T',  1._wp )
855
856         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
857         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
858#if defined key_agrif
859         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
860#endif
861         
862         ! Sea Surface Height at u-,v-points (vvl case only)
863         IF( .NOT.ln_linssh ) THEN                               
864            DO jj = 2, jpjm1
865               DO ji = 2, jpim1      ! NO Vector Opt.
866                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
867                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
868                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
869                  zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
870                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
871                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
872               END DO
873            END DO
874            CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
875         ENDIF   
876         !                                 
877         ! Half-step back interpolation of SSH for surface pressure computation:
878         !----------------------------------------------------------------------
879         IF ((jn==1).AND.ll_init) THEN
880           za0=1._wp                        ! Forward-backward
881           za1=0._wp                           
882           za2=0._wp
883           za3=0._wp
884         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
885           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
886           za1=-0.1666666666666_wp          ! za1 = gam
887           za2= 0.0833333333333_wp          ! za2 = eps
888           za3= 0._wp             
889         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
890            IF (rn_bt_alpha==0._wp) THEN
891               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps
892               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps
893               za2=0.088_wp                 ! za2 = gam
894               za3=0.013_wp                 ! za3 = eps
895            ELSE
896               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
897               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
898               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
899               za1 = 1._wp - za0 - zgamma - zepsilon
900               za2 = zgamma
901               za3 = zepsilon
902            ENDIF
903         ENDIF
904         !
905         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
906            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
907         
908         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters
909           DO jj = 2, jpjm1
910              DO ji = 2, jpim1 
911                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
912                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            &
913                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
914                     &                                                             > rn_wdmin1 + rn_wdmin2
915                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
916                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
917                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
918   
919                IF(ll_tmp1) THEN
920                  zcpx(ji,jj) = 1.0_wp
921                ELSE IF(ll_tmp2) THEN
922                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here
923                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
924                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) )
925                ELSE
926                  zcpx(ji,jj) = 0._wp
927                ENDIF
928                !
929                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
930                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            &
931                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
932                     &                                                             > rn_wdmin1 + rn_wdmin2
933                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
934                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
935                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
936   
937                IF(ll_tmp1) THEN
938                  zcpy(ji,jj) = 1.0_wp
939                ELSEIF(ll_tmp2) THEN
940                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
941                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
942                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) )
943                ELSE
944                  zcpy(ji,jj) = 0._wp
945                ENDIF
946              END DO
947           END DO
948         ENDIF
949         !
950         ! Compute associated depths at U and V points:
951         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form
952            !                                       
953            DO jj = 2, jpjm1                           
954               DO ji = 2, jpim1
955                  zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
956                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
957                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
958                  zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
959                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
960                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
961                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
962                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
963               END DO
964            END DO
965            !
966         ENDIF
967         !
968         ! Add Coriolis trend:
969         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
970         ! at each time step. We however keep them constant here for optimization.
971         ! Recall that zwx and zwy arrays hold fluxes at this stage:
972         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
973         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
974         !
975         SELECT CASE( nvor_scheme )
976         CASE( np_ENT )             ! energy conserving scheme (t-point)
977         DO jj = 2, jpjm1
978            DO ji = 2, jpim1   ! vector opt.
979
980               z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
981               z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
982           
983               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   &
984                  &               * (  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) )   &
985                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   )
986                  !
987               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
988                  &               * (  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) )   & 
989                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   ) 
990            END DO 
991         END DO 
992         !         
993         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point)
994            DO jj = 2, jpjm1
995               DO ji = fs_2, fs_jpim1   ! vector opt.
996                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
997                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
998                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
999                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
1000                  zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1001                  zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1002               END DO
1003            END DO
1004            !
1005         CASE( np_ENS )             ! enstrophy conserving scheme (f-point)
1006            DO jj = 2, jpjm1
1007               DO ji = fs_2, fs_jpim1   ! vector opt.
1008                  zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
1009                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
1010                  zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
1011                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
1012                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1013                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1014               END DO
1015            END DO
1016            !
1017         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f)
1018            DO jj = 2, jpjm1
1019               DO ji = fs_2, fs_jpim1   ! vector opt.
1020                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  &
1021                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  &
1022                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  & 
1023                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  )
1024                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  & 
1025                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  &
1026                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  & 
1027                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  )
1028               END DO
1029            END DO
1030            !
1031         END SELECT
1032         !
1033         ! Add tidal astronomical forcing if defined
1034         IF ( ln_tide .AND. ln_tide_pot ) THEN
1035            DO jj = 2, jpjm1
1036               DO ji = fs_2, fs_jpim1   ! vector opt.
1037                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1038                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1039                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1040                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1041               END DO
1042            END DO
1043         ENDIF
1044         !
1045         ! Add bottom stresses:
1046!jth do implicitly instead
1047         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1048            DO jj = 2, jpjm1
1049               DO ji = fs_2, fs_jpim1   ! vector opt.
1050                  zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1051                  zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1052               END DO
1053            END DO
1054         ENDIF 
1055         !
1056         ! Surface pressure trend:
1057         IF( ln_wd_il ) THEN
1058           DO jj = 2, jpjm1
1059              DO ji = 2, jpim1 
1060                 ! Add surface pressure gradient
1061                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1062                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1063                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 
1064                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1065              END DO
1066           END DO
1067         ELSE
1068           DO jj = 2, jpjm1
1069              DO ji = fs_2, fs_jpim1   ! vector opt.
1070                 ! Add surface pressure gradient
1071                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1072                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1073                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1074                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1075              END DO
1076           END DO
1077         END IF
1078
1079         !
1080         ! Set next velocities:
1081         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
1082            DO jj = 2, jpjm1
1083               DO ji = fs_2, fs_jpim1   ! vector opt.
1084                  ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
1085                            &     + rdtbt * (                      zwx(ji,jj)   &
1086                            &                                 + zu_trd(ji,jj)   &
1087                            &                                 + zu_frc(ji,jj) ) & 
1088                            &   ) * ssumask(ji,jj)
1089
1090                  va_e(ji,jj) = (                                 vn_e(ji,jj)   &
1091                            &     + rdtbt * (                      zwy(ji,jj)   &
1092                            &                                 + zv_trd(ji,jj)   &
1093                            &                                 + zv_frc(ji,jj) ) &
1094                            &   ) * ssvmask(ji,jj)
1095 
1096               END DO
1097            END DO
1098            !
1099         ELSE                           !* Flux form
1100            DO jj = 2, jpjm1
1101               DO ji = fs_2, fs_jpim1   ! vector opt.
1102
1103                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1104                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1105
1106                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1107                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1108
1109                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
1110                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
1111                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
1112                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
1113                            &   ) * zhura
1114
1115                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
1116                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
1117                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
1118                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
1119                            &   ) * zhvra
1120               END DO
1121            END DO
1122         ENDIF
1123!jth implicit bottom friction:
1124         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1125            DO jj = 2, jpjm1
1126               DO ji = fs_2, fs_jpim1   ! vector opt.
1127                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1128                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1129               END DO
1130            END DO
1131         ENDIF
1132
1133         
1134         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
1135            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1136            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1137            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1138            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1139            !
1140         ENDIF
1141         !                                             !* domain lateral boundary
1142         CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1143         !
1144         !                                                 ! open boundaries
1145         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1146#if defined key_agrif                                                           
1147         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
1148#endif
1149         !                                             !* Swap
1150         !                                             !  ----
1151         ubb_e  (:,:) = ub_e  (:,:)
1152         ub_e   (:,:) = un_e  (:,:)
1153         un_e   (:,:) = ua_e  (:,:)
1154         !
1155         vbb_e  (:,:) = vb_e  (:,:)
1156         vb_e   (:,:) = vn_e  (:,:)
1157         vn_e   (:,:) = va_e  (:,:)
1158         !
1159         sshbb_e(:,:) = sshb_e(:,:)
1160         sshb_e (:,:) = sshn_e(:,:)
1161         sshn_e (:,:) = ssha_e(:,:)
1162
1163         !                                             !* Sum over whole bt loop
1164         !                                             !  ----------------------
1165         za1 = wgtbtp1(jn)                                   
1166         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
1167            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
1168            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
1169         ELSE                                       ! Sum transports
1170            IF ( .NOT.ln_wd_dl ) THEN 
1171               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
1172               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
1173            ELSE
1174               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
1175               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
1176            END IF
1177         ENDIF
1178         !                                          ! Sum sea level
1179         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1180
1181         !                                                 ! ==================== !
1182      END DO                                               !        end loop      !
1183      !                                                    ! ==================== !
1184      ! -----------------------------------------------------------------------------
1185      ! Phase 3. update the general trend with the barotropic trend
1186      ! -----------------------------------------------------------------------------
1187      !
1188      ! Set advection velocity correction:
1189      IF (ln_bt_fw) THEN
1190         zwx(:,:) = un_adv(:,:)
1191         zwy(:,:) = vn_adv(:,:)
1192         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1193            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1194            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1195            !
1196            ! Update corrective fluxes for next time step:
1197            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1198            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1199         ELSE
1200            un_bf(:,:) = 0._wp
1201            vn_bf(:,:) = 0._wp 
1202         END IF         
1203         ! Save integrated transport for next computation
1204         ub2_b(:,:) = zwx(:,:)
1205         vb2_b(:,:) = zwy(:,:)
1206      ENDIF
1207
1208
1209      !
1210      ! Update barotropic trend:
1211      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1212         DO jk=1,jpkm1
1213            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1214            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1215         END DO
1216      ELSE
1217         ! At this stage, ssha has been corrected: compute new depths at velocity points
1218         DO jj = 1, jpjm1
1219            DO ji = 1, jpim1      ! NO Vector Opt.
1220               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
1221                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      &
1222                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1223               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
1224                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      &
1225                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1226            END DO
1227         END DO
1228         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1229         !
1230         DO jk=1,jpkm1
1231            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1232            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1233         END DO
1234         ! Save barotropic velocities not transport:
1235         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1236         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1237      ENDIF
1238
1239
1240      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
1241      DO jk = 1, jpkm1
1242         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)
1243         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1244      END DO
1245
1246      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1247         DO jk = 1, jpkm1
1248            un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1249                       & + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk) 
1250            vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) & 
1251                       & + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk) 
1252         END DO
1253      END IF
1254
1255     
1256      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu_n(:,:) )    ! barotropic i-current
1257      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv_n(:,:) )    ! barotropic i-current
1258      !
1259#if defined key_agrif
1260      ! Save time integrated fluxes during child grid integration
1261      ! (used to update coarse grid transports at next time step)
1262      !
1263      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1264         IF( Agrif_NbStepint() == 0 ) THEN
1265            ub2_i_b(:,:) = 0._wp
1266            vb2_i_b(:,:) = 0._wp
1267         END IF
1268         !
1269         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1270         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1271         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1272      ENDIF
1273#endif     
1274      !                                   !* write time-spliting arrays in the restart
1275      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
1276      !
1277      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
1278      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1279      !
1280      IF( ln_diatmb ) THEN
1281         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
1282         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
1283      ENDIF
1284      !
1285   END SUBROUTINE dyn_spg_ts
1286
1287
1288   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1289      !!---------------------------------------------------------------------
1290      !!                   ***  ROUTINE ts_wgt  ***
1291      !!
1292      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1293      !!----------------------------------------------------------------------
1294      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
1295      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
1296      INTEGER, INTENT(inout) :: jpit      ! cycle length   
1297      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
1298                                                         zwgt2    ! Secondary weights
1299     
1300      INTEGER ::  jic, jn, ji                      ! temporary integers
1301      REAL(wp) :: za1, za2
1302      !!----------------------------------------------------------------------
1303
1304      zwgt1(:) = 0._wp
1305      zwgt2(:) = 0._wp
1306
1307      ! Set time index when averaged value is requested
1308      IF (ll_fw) THEN
1309         jic = nn_baro
1310      ELSE
1311         jic = 2 * nn_baro
1312      ENDIF
1313
1314      ! Set primary weights:
1315      IF (ll_av) THEN
1316           ! Define simple boxcar window for primary weights
1317           ! (width = nn_baro, centered around jic)     
1318         SELECT CASE ( nn_bt_flt )
1319              CASE( 0 )  ! No averaging
1320                 zwgt1(jic) = 1._wp
1321                 jpit = jic
1322
1323              CASE( 1 )  ! Boxcar, width = nn_baro
1324                 DO jn = 1, 3*nn_baro
1325                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1326                    IF (za1 < 0.5_wp) THEN
1327                      zwgt1(jn) = 1._wp
1328                      jpit = jn
1329                    ENDIF
1330                 ENDDO
1331
1332              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
1333                 DO jn = 1, 3*nn_baro
1334                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1335                    IF (za1 < 1._wp) THEN
1336                      zwgt1(jn) = 1._wp
1337                      jpit = jn
1338                    ENDIF
1339                 ENDDO
1340              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1341         END SELECT
1342
1343      ELSE ! No time averaging
1344         zwgt1(jic) = 1._wp
1345         jpit = jic
1346      ENDIF
1347   
1348      ! Set secondary weights
1349      DO jn = 1, jpit
1350        DO ji = jn, jpit
1351             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1352        END DO
1353      END DO
1354
1355      ! Normalize weigths:
1356      za1 = 1._wp / SUM(zwgt1(1:jpit))
1357      za2 = 1._wp / SUM(zwgt2(1:jpit))
1358      DO jn = 1, jpit
1359        zwgt1(jn) = zwgt1(jn) * za1
1360        zwgt2(jn) = zwgt2(jn) * za2
1361      END DO
1362      !
1363   END SUBROUTINE ts_wgt
1364
1365
1366   SUBROUTINE ts_rst( kt, cdrw )
1367      !!---------------------------------------------------------------------
1368      !!                   ***  ROUTINE ts_rst  ***
1369      !!
1370      !! ** Purpose : Read or write time-splitting arrays in restart file
1371      !!----------------------------------------------------------------------
1372      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
1373      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1374      !!----------------------------------------------------------------------
1375      !
1376      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1377         !                                   ! ---------------
1378         IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN    !* Read the restart file
1379            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
1380            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
1381            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
1382            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
1383            IF( .NOT.ln_bt_av ) THEN
1384               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
1385               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
1386               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
1387               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
1388               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
1389               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
1390            ENDIF
1391#if defined key_agrif
1392            ! Read time integrated fluxes
1393            IF ( .NOT.Agrif_Root() ) THEN
1394               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
1395               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
1396            ENDIF
1397#endif
1398         ELSE                                   !* Start from rest
1399            IF(lwp) WRITE(numout,*)
1400            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
1401            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
1402            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
1403            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
1404#if defined key_agrif
1405            IF ( .NOT.Agrif_Root() ) THEN
1406               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
1407            ENDIF
1408#endif
1409         ENDIF
1410         !
1411      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1412         !                                   ! -------------------
1413         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1414         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1415         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
1416         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
1417         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
1418         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
1419         !
1420         IF (.NOT.ln_bt_av) THEN
1421            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
1422            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
1423            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
1424            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
1425            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
1426            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
1427         ENDIF
1428#if defined key_agrif
1429         ! Save time integrated fluxes
1430         IF ( .NOT.Agrif_Root() ) THEN
1431            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
1432            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
1433         ENDIF
1434#endif
1435         IF( lwxios ) CALL iom_swap(      cxios_context          )
1436      ENDIF
1437      !
1438   END SUBROUTINE ts_rst
1439
1440
1441   SUBROUTINE dyn_spg_ts_init
1442      !!---------------------------------------------------------------------
1443      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1444      !!
1445      !! ** Purpose : Set time splitting options
1446      !!----------------------------------------------------------------------
1447      INTEGER  ::   ji ,jj              ! dummy loop indices
1448      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1449      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
1450      !!----------------------------------------------------------------------
1451      !
1452      ! Max courant number for ext. grav. waves
1453      !
1454      DO jj = 1, jpj
1455         DO ji =1, jpi
1456            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1457            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1458            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1459         END DO
1460      END DO
1461      !
1462      zcmax = MAXVAL( zcu(:,:) )
1463      CALL mpp_max( 'dynspg_ts', zcmax )
1464
1465      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1466      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1467     
1468      rdtbt = rdt / REAL( nn_baro , wp )
1469      zcmax = zcmax * rdtbt
1470      ! Print results
1471      IF(lwp) WRITE(numout,*)
1472      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1473      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1474      IF( ln_bt_auto ) THEN
1475         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
1476         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1477      ELSE
1478         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
1479      ENDIF
1480
1481      IF(ln_bt_av) THEN
1482         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
1483      ELSE
1484         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1485      ENDIF
1486      !
1487      !
1488      IF(ln_bt_fw) THEN
1489         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1490      ELSE
1491         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1492      ENDIF
1493      !
1494#if defined key_agrif
1495      ! Restrict the use of Agrif to the forward case only
1496!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1497#endif
1498      !
1499      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1500      SELECT CASE ( nn_bt_flt )
1501         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1502         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1503         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1504         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1505      END SELECT
1506      !
1507      IF(lwp) WRITE(numout,*) ' '
1508      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1509      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1510      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1511      !
1512      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1513      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1514         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1515      ENDIF
1516      !
1517      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1518         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1519      ENDIF
1520      IF( zcmax>0.9_wp ) THEN
1521         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1522      ENDIF
1523      !
1524      !                             ! Allocate time-splitting arrays
1525      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1526      !
1527      !                             ! read restart when needed
1528      CALL ts_rst( nit000, 'READ' )
1529      !
1530      IF( lwxios ) THEN
1531! define variables in restart file when writing with XIOS
1532         CALL iom_set_rstw_var_active('ub2_b')
1533         CALL iom_set_rstw_var_active('vb2_b')
1534         CALL iom_set_rstw_var_active('un_bf')
1535         CALL iom_set_rstw_var_active('vn_bf')
1536         !
1537         IF (.NOT.ln_bt_av) THEN
1538            CALL iom_set_rstw_var_active('sshbb_e')
1539            CALL iom_set_rstw_var_active('ubb_e')
1540            CALL iom_set_rstw_var_active('vbb_e')
1541            CALL iom_set_rstw_var_active('sshb_e')
1542            CALL iom_set_rstw_var_active('ub_e')
1543            CALL iom_set_rstw_var_active('vb_e')
1544         ENDIF
1545#if defined key_agrif
1546         ! Save time integrated fluxes
1547         IF ( .NOT.Agrif_Root() ) THEN
1548            CALL iom_set_rstw_var_active('ub2_i_b')
1549            CALL iom_set_rstw_var_active('vb2_i_b')
1550         ENDIF
1551#endif
1552      ENDIF
1553      !
1554   END SUBROUTINE dyn_spg_ts_init
1555
1556   !!======================================================================
1557END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.