source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 13289

Last change on this file since 13289 was 13289, checked in by smasson, 3 months ago

trunk: dynspg_ts: must call lnc_lnk before ln_bdy staff

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