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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9506

Last change on this file since 9506 was 9506, checked in by clem, 4 years ago

for some mpp decomposition (typically 8x4) agrif failed at running because some barotropic fields needed to be initialized (ub2_b etc) for the 1st interpolation and update. These fields are now set to 0 in dynspg_ts

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