source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6986

Last change on this file since 6986 was 6986, checked in by acc, 5 years ago

Branch dev_r6393_NOC_WAD. Introduced some WAD test cases, tidied and corrected WAD code

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