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_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

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

Last change on this file since 8842 was 8842, checked in by jchanut, 6 years ago

Enhance4-freesurface. step 3: diffusive barotropic time stepping - #1959

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