source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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