source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_ts.F90 @ 9863

Last change on this file since 9863 was 9863, checked in by gm, 3 years ago

#1911 (ENHANCE-04): simplified implementation of the Euler stepping at nit000

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