New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in branches/2017/dev_r8600_xios_read/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_r8600_xios_read/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8612

Last change on this file since 8612 was 8612, checked in by andmirek, 6 years ago

#1953 read single file restart with XIOS

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