source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8576

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

#1883 (HPC-09) - correction of a bug in dynspg_ts with top drag

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