source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 4 years ago

update trunk with OpenMP parallelization

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