source: branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8888

Last change on this file since 8888 was 8888, checked in by jchanut, 3 years ago

Enhance4-freesurface. Minor corrections after testing, update namelists - #1959

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