source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7831

Last change on this file since 7831 was 7831, checked in by jamesharle, 3 years ago

Remove key_tide from dynspg_ts.F90 (ticket #1874)

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