source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 12732

Last change on this file since 12732 was 12732, checked in by techene, 6 months ago

some cleaning and proper module/routine name, mini bug introduced and corrected in sbcice_cice

  • Property svn:keywords set to Id
File size: 78.5 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         !                             ! Sum over sub-time-steps to compute advective velocities
525         za2 = wgtbtp2(jn)             ! zhU, zhV hold fluxes extrapolated at jn+0.5
526         un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
527         vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
528         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
529         IF ( ln_wd_dl_bc ) THEN
530            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
531            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row
532         END IF
533         !
534         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
535         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
536#if defined key_agrif
537         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
538#endif
539         
540         ! Sea Surface Height at u-,v-points (vvl case only)
541         IF( .NOT.ln_linssh ) THEN                               
542            DO_2D_00_00
543               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
544                  &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
545                  &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
546               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
547                  &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
548                  &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
549            END_2D
550         ENDIF   
551         !         
552         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
553         !--            m+1/2           m+1              m               m-1              m-2     --!
554         !--        ssh'    =  za0 * ssh     +  za1 * ssh   +  za2 * ssh      +  za3 * ssh        --!
555         !------------------------------------------------------------------------------------------!
556         CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 )   ! coeficients of the interpolation
557         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
558            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
559         !
560         !                             ! Surface pressure gradient
561         zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
562         DO_2D_00_00
563            zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
564            zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
565         END_2D
566         IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient
567            CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters
568            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
569            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
570         ENDIF
571         !
572         ! Add Coriolis trend:
573         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
574         ! at each time step. We however keep them constant here for optimization.
575         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
576         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
577         !
578         ! Add tidal astronomical forcing if defined
579         IF ( ln_tide .AND. ln_tide_pot ) THEN
580            DO_2D_00_00
581               zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
582               zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
583            END_2D
584         ENDIF
585         !
586         ! Add bottom stresses:
587!jth do implicitly instead
588         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
589            DO_2D_00_00
590               zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
591               zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
592            END_2D
593         ENDIF
594         !
595         ! Set next velocities:
596         !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn)
597         !--                              VECTOR FORM
598         !--   m+1                 m               /                                                       m+1/2           \    --!
599         !--  u     =             u   + delta_t' * \         (1-r)*g * grad_x( ssh') -         f * k vect u      +     frc /    --!
600         !--                                                                                                                    --!
601         !--                             FLUX FORM                                                                              --!
602         !--  m+1   __1__  /  m    m               /  m+1/2                             m+1/2              m+1/2    n      \ \  --!
603         !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --!
604         !--         h     \                                                                                                 /  --!
605         !------------------------------------------------------------------------------------------------------------------------!
606         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
607            DO_2D_00_00
608               ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
609                         &     + rDt_e * (                   zu_spg(ji,jj)   &
610                         &                                 + zu_trd(ji,jj)   &
611                         &                                 + zu_frc(ji,jj) ) & 
612                         &   ) * ssumask(ji,jj)
613
614               va_e(ji,jj) = (                                 vn_e(ji,jj)   &
615                         &     + rDt_e * (                   zv_spg(ji,jj)   &
616                         &                                 + zv_trd(ji,jj)   &
617                         &                                 + zv_frc(ji,jj) ) &
618                         &   ) * ssvmask(ji,jj)
619            END_2D
620            !
621         ELSE                           !* Flux form
622            DO_2D_00_00
623               !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
624               !                    ! backward interpolated depth used in spg terms at jn+1/2
625               zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
626                    &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
627               zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
628                    &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
629               !                    ! inverse depth at jn+1
630               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
631               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
632               !
633               ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
634                    &            + rDt_e * (  zhu_bck        * zu_spg (ji,jj)  &   !
635                    &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
636                    &                       +  hu(ji,jj,Kmm) * zu_frc (ji,jj)  )   ) * z1_hu
637               !
638               va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      &
639                    &            + rDt_e * (  zhv_bck        * zv_spg (ji,jj)  &   !
640                    &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   !
641                    &                       +  hv(ji,jj,Kmm) * zv_frc (ji,jj)  )   ) * z1_hv
642            END_2D
643         ENDIF
644!jth implicit bottom friction:
645         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
646            DO_2D_00_00
647                  ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
648                  va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
649            END_2D
650         ENDIF
651       
652         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
653            hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
654            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
655            hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
656            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
657            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
658                 &                         , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
659                 &                         , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
660         ELSE
661            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
662         ENDIF
663         !
664         !
665         !                                                 ! open boundaries
666         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
667#if defined key_agrif                                                           
668         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
669#endif
670         !                                             !* Swap
671         !                                             !  ----
672         ubb_e  (:,:) = ub_e  (:,:)
673         ub_e   (:,:) = un_e  (:,:)
674         un_e   (:,:) = ua_e  (:,:)
675         !
676         vbb_e  (:,:) = vb_e  (:,:)
677         vb_e   (:,:) = vn_e  (:,:)
678         vn_e   (:,:) = va_e  (:,:)
679         !
680         sshbb_e(:,:) = sshb_e(:,:)
681         sshb_e (:,:) = sshn_e(:,:)
682         sshn_e (:,:) = ssha_e(:,:)
683
684         !                                             !* Sum over whole bt loop
685         !                                             !  ----------------------
686         za1 = wgtbtp1(jn)                                   
687         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
688            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
689            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
690         ELSE                                       ! Sum transports
691            IF ( .NOT.ln_wd_dl ) THEN 
692               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
693               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
694            ELSE
695               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
696               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
697            END IF
698         ENDIF
699         !                                          ! Sum sea level
700         pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
701
702         !                                                 ! ==================== !
703      END DO                                               !        end loop      !
704      !                                                    ! ==================== !
705      ! -----------------------------------------------------------------------------
706      ! Phase 3. update the general trend with the barotropic trend
707      ! -----------------------------------------------------------------------------
708      !
709      ! Set advection velocity correction:
710      IF (ln_bt_fw) THEN
711         IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
712            DO_2D_11_11
713               zun_save = un_adv(ji,jj)
714               zvn_save = vn_adv(ji,jj)
715               !                          ! apply the previously computed correction
716               un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
717               vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
718               !                          ! Update corrective fluxes for next time step
719               un_bf(ji,jj)  = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
720               vn_bf(ji,jj)  = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
721               !                          ! Save integrated transport for next computation
722               ub2_b(ji,jj) = zun_save
723               vb2_b(ji,jj) = zvn_save
724            END_2D
725         ELSE
726            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero
727            vn_bf(:,:) = 0._wp
728            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
729            vb2_b(:,:) = vn_adv(:,:)
730         END IF
731      ENDIF
732
733
734      !
735      ! Update barotropic trend:
736      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
737         DO jk=1,jpkm1
738            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
739            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
740         END DO
741      ELSE
742         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
743         DO_2D_10_10
744            zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
745               &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      &
746               &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
747            zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
748               &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      &
749               &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
750         END_2D
751         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
752         !
753         DO jk=1,jpkm1
754            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
755            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
756         END DO
757         ! Save barotropic velocities not transport:
758         puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
759         pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
760      ENDIF
761
762
763      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
764      DO jk = 1, jpkm1
765         puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
766         pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)
767      END DO
768
769      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
770         DO jk = 1, jpkm1
771            puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
772                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
773            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
774                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 
775         END DO
776      END IF
777
778     
779      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
780      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
781      !
782#if defined key_agrif
783      ! Save time integrated fluxes during child grid integration
784      ! (used to update coarse grid transports at next time step)
785      !
786      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
787         IF( Agrif_NbStepint() == 0 ) THEN
788            ub2_i_b(:,:) = 0._wp
789            vb2_i_b(:,:) = 0._wp
790         END IF
791         !
792         za1 = 1._wp / REAL(Agrif_rhot(), wp)
793         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
794         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
795      ENDIF
796#endif     
797      !                                   !* write time-spliting arrays in the restart
798      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
799      !
800      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
801      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
802      !
803      CALL iom_put( "baro_u" , puu_b(:,:,Kmm) )  ! Barotropic  U Velocity
804      CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) )  ! Barotropic  V Velocity
805      !
806   END SUBROUTINE dyn_spg_ts
807
808
809   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
810      !!---------------------------------------------------------------------
811      !!                   ***  ROUTINE ts_wgt  ***
812      !!
813      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
814      !!----------------------------------------------------------------------
815      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
816      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
817      INTEGER, INTENT(inout) :: jpit      ! cycle length   
818      REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights
819                                                         zwgt2    ! Secondary weights
820     
821      INTEGER ::  jic, jn, ji                      ! temporary integers
822      REAL(wp) :: za1, za2
823      !!----------------------------------------------------------------------
824
825      zwgt1(:) = 0._wp
826      zwgt2(:) = 0._wp
827
828      ! Set time index when averaged value is requested
829      IF (ll_fw) THEN
830         jic = nn_e
831      ELSE
832         jic = 2 * nn_e
833      ENDIF
834
835      ! Set primary weights:
836      IF (ll_av) THEN
837           ! Define simple boxcar window for primary weights
838           ! (width = nn_e, centered around jic)     
839         SELECT CASE ( nn_bt_flt )
840              CASE( 0 )  ! No averaging
841                 zwgt1(jic) = 1._wp
842                 jpit = jic
843
844              CASE( 1 )  ! Boxcar, width = nn_e
845                 DO jn = 1, 3*nn_e
846                    za1 = ABS(float(jn-jic))/float(nn_e) 
847                    IF (za1 < 0.5_wp) THEN
848                      zwgt1(jn) = 1._wp
849                      jpit = jn
850                    ENDIF
851                 ENDDO
852
853              CASE( 2 )  ! Boxcar, width = 2 * nn_e
854                 DO jn = 1, 3*nn_e
855                    za1 = ABS(float(jn-jic))/float(nn_e) 
856                    IF (za1 < 1._wp) THEN
857                      zwgt1(jn) = 1._wp
858                      jpit = jn
859                    ENDIF
860                 ENDDO
861              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
862         END SELECT
863
864      ELSE ! No time averaging
865         zwgt1(jic) = 1._wp
866         jpit = jic
867      ENDIF
868   
869      ! Set secondary weights
870      DO jn = 1, jpit
871        DO ji = jn, jpit
872             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
873        END DO
874      END DO
875
876      ! Normalize weigths:
877      za1 = 1._wp / SUM(zwgt1(1:jpit))
878      za2 = 1._wp / SUM(zwgt2(1:jpit))
879      DO jn = 1, jpit
880        zwgt1(jn) = zwgt1(jn) * za1
881        zwgt2(jn) = zwgt2(jn) * za2
882      END DO
883      !
884   END SUBROUTINE ts_wgt
885
886
887   SUBROUTINE ts_rst( kt, cdrw )
888      !!---------------------------------------------------------------------
889      !!                   ***  ROUTINE ts_rst  ***
890      !!
891      !! ** Purpose : Read or write time-splitting arrays in restart file
892      !!----------------------------------------------------------------------
893      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
894      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
895      !!----------------------------------------------------------------------
896      !
897      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
898         !                                   ! ---------------
899         IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN    !* Read the restart file
900            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
901            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
902            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
903            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
904            IF( .NOT.ln_bt_av ) THEN
905               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
906               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
907               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
908               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
909               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
910               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
911            ENDIF
912#if defined key_agrif
913            ! Read time integrated fluxes
914            IF ( .NOT.Agrif_Root() ) THEN
915               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
916               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
917            ENDIF
918#endif
919         ELSE                                   !* Start from rest
920            IF(lwp) WRITE(numout,*)
921            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
922            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
923            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
924            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
925#if defined key_agrif
926            IF ( .NOT.Agrif_Root() ) THEN
927               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
928            ENDIF
929#endif
930         ENDIF
931         !
932      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
933         !                                   ! -------------------
934         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
935         IF( lwxios ) CALL iom_swap(      cwxios_context          )
936         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
937         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
938         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
939         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
940         !
941         IF (.NOT.ln_bt_av) THEN
942            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
943            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
944            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
945            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
946            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
947            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
948         ENDIF
949#if defined key_agrif
950         ! Save time integrated fluxes
951         IF ( .NOT.Agrif_Root() ) THEN
952            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
953            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
954         ENDIF
955#endif
956         IF( lwxios ) CALL iom_swap(      cxios_context          )
957      ENDIF
958      !
959   END SUBROUTINE ts_rst
960
961
962   SUBROUTINE dyn_spg_ts_init
963      !!---------------------------------------------------------------------
964      !!                   ***  ROUTINE dyn_spg_ts_init  ***
965      !!
966      !! ** Purpose : Set time splitting options
967      !!----------------------------------------------------------------------
968      INTEGER  ::   ji ,jj              ! dummy loop indices
969      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
970      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
971      !!----------------------------------------------------------------------
972      !
973      ! Max courant number for ext. grav. waves
974      !
975      DO_2D_11_11
976         zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
977         zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
978         zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
979      END_2D
980      !
981      zcmax = MAXVAL( zcu(:,:) )
982      CALL mpp_max( 'dynspg_ts', zcmax )
983
984      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
985      IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
986     
987      rDt_e = rn_Dt / REAL( nn_e , wp )
988      zcmax = zcmax * rDt_e
989      ! Print results
990      IF(lwp) WRITE(numout,*)
991      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
992      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
993      IF( ln_bt_auto ) THEN
994         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_e '
995         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
996      ELSE
997         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e
998      ENDIF
999
1000      IF(ln_bt_av) THEN
1001         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on '
1002      ELSE
1003         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1004      ENDIF
1005      !
1006      !
1007      IF(ln_bt_fw) THEN
1008         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1009      ELSE
1010         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1011      ENDIF
1012      !
1013#if defined key_agrif
1014      ! Restrict the use of Agrif to the forward case only
1015!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1016#endif
1017      !
1018      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1019      SELECT CASE ( nn_bt_flt )
1020         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1021         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
1022         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
1023         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1024      END SELECT
1025      !
1026      IF(lwp) WRITE(numout,*) ' '
1027      IF(lwp) WRITE(numout,*) '     nn_e = ', nn_e
1028      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rDt_e
1029      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1030      !
1031      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1032      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1033         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1034      ENDIF
1035      !
1036      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1037         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1038      ENDIF
1039      IF( zcmax>0.9_wp ) THEN
1040         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )         
1041      ENDIF
1042      !
1043      !                             ! Allocate time-splitting arrays
1044      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1045      !
1046      !                             ! read restart when needed
1047      CALL ts_rst( nit000, 'READ' )
1048      !
1049      IF( lwxios ) THEN
1050! define variables in restart file when writing with XIOS
1051         CALL iom_set_rstw_var_active('ub2_b')
1052         CALL iom_set_rstw_var_active('vb2_b')
1053         CALL iom_set_rstw_var_active('un_bf')
1054         CALL iom_set_rstw_var_active('vn_bf')
1055         !
1056         IF (.NOT.ln_bt_av) THEN
1057            CALL iom_set_rstw_var_active('sshbb_e')
1058            CALL iom_set_rstw_var_active('ubb_e')
1059            CALL iom_set_rstw_var_active('vbb_e')
1060            CALL iom_set_rstw_var_active('sshb_e')
1061            CALL iom_set_rstw_var_active('ub_e')
1062            CALL iom_set_rstw_var_active('vb_e')
1063         ENDIF
1064#if defined key_agrif
1065         ! Save time integrated fluxes
1066         IF ( .NOT.Agrif_Root() ) THEN
1067            CALL iom_set_rstw_var_active('ub2_i_b')
1068            CALL iom_set_rstw_var_active('vb2_i_b')
1069         ENDIF
1070#endif
1071      ENDIF
1072      !
1073   END SUBROUTINE dyn_spg_ts_init
1074
1075   
1076   SUBROUTINE dyn_cor_2D_init( Kmm )
1077      !!---------------------------------------------------------------------
1078      !!                   ***  ROUTINE dyn_cor_2D_init  ***
1079      !!
1080      !! ** Purpose : Set time splitting options
1081      !! Set arrays to remove/compute coriolis trend.
1082      !! Do it once during initialization if volume is fixed, else at each long time step.
1083      !! Note that these arrays are also used during barotropic loop. These are however frozen
1084      !! although they should be updated in the variable volume case. Not a big approximation.
1085      !! To remove this approximation, copy lines below inside barotropic loop
1086      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1087      !!
1088      !! Compute zwz = f / ( height of the water colomn )
1089      !!----------------------------------------------------------------------
1090      INTEGER,  INTENT(in)         ::  Kmm  ! Time index
1091      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
1092      REAL(wp) ::   z1_ht
1093      REAL(wp), DIMENSION(jpi,jpj) :: zhf
1094      !!----------------------------------------------------------------------
1095      !
1096      SELECT CASE( nvor_scheme )
1097      CASE( np_EEN )                != EEN scheme using e3f energy & enstrophy scheme
1098         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
1099         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
1100            DO_2D_10_10
1101               zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
1102                    &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp 
1103               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1104            END_2D
1105         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
1106            DO_2D_10_10
1107               zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
1108                    &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
1109                    &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
1110                    &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
1111               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1112            END_2D
1113         END SELECT
1114         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1115         !
1116         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1117         DO_2D_01_01
1118            ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1119            ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1120            ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1121            ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1122         END_2D
1123         !
1124      CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme
1125         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1126         DO_2D_01_01
1127            z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1128            ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
1129            ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
1130            ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1131            ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
1132         END_2D
1133         !
1134      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
1135         !
1136         zwz(:,:) = 0._wp
1137         zhf(:,:) = 0._wp
1138         
1139         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1140!!gm    A priori a better value should be something like :
1141!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1142!!gm                     divided by the sum of the corresponding mask
1143!!gm
1144!!           
1145         IF( .NOT.ln_sco ) THEN
1146 
1147   !!gm  agree the JC comment  : this should be done in a much clear way
1148 
1149   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1150   !     Set it to zero for the time being
1151   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
1152   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
1153   !              ENDIF
1154   !              zhf(:,:) = gdepw_0(:,:,jk+1)
1155            !
1156         ELSE
1157            !
1158            !zhf(:,:) = hbatf(:,:)
1159            DO_2D_10_10
1160               zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
1161                    &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
1162                    &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
1163                    &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
1164            END_2D
1165         ENDIF
1166         !
1167         DO jj = 1, jpjm1
1168            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1169         END DO
1170         !
1171         DO jk = 1, jpkm1
1172            DO jj = 1, jpjm1
1173               zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1174            END DO
1175         END DO
1176         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1177         ! JC: TBC. hf should be greater than 0
1178         DO_2D_11_11
1179            IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
1180         END_2D
1181         zwz(:,:) = ff_f(:,:) * zwz(:,:)
1182      END SELECT
1183     
1184   END SUBROUTINE dyn_cor_2d_init
1185
1186
1187
1188   SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
1189      !!---------------------------------------------------------------------
1190      !!                   ***  ROUTINE dyn_cor_2d  ***
1191      !!
1192      !! ** Purpose : Compute u and v coriolis trends
1193      !!----------------------------------------------------------------------
1194      INTEGER  ::   ji ,jj                             ! dummy loop indices
1195      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
1196      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pht, phu, phv, punb, pvnb, zhU, zhV
1197      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
1198      !!----------------------------------------------------------------------
1199      SELECT CASE( nvor_scheme )
1200      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
1201         DO_2D_00_00
1202            z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1203            z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1204            zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
1205               &               * (  e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
1206               &                  + e1e2t(ji  ,jj)*pht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
1207               !
1208            zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1209               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
1210               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
1211         END_2D
1212         !         
1213      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
1214         DO_2D_00_00
1215            zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1216            zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1217            zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1218            zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1219            ! energy conserving formulation for planetary vorticity term
1220            zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1221            zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1222         END_2D
1223         !
1224      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
1225         DO_2D_00_00
1226            zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
1227              &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1228            zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
1229              &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1230            zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1231            zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1232         END_2D
1233         !
1234      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
1235         DO_2D_00_00
1236            zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
1237             &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
1238             &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
1239             &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1240            zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1241             &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
1242             &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
1243             &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
1244         END_2D
1245         !
1246      END SELECT
1247      !
1248   END SUBROUTINE dyn_cor_2D
1249
1250
1251   SUBROUTINE wad_tmsk( pssh, ptmsk )
1252      !!----------------------------------------------------------------------
1253      !!                  ***  ROUTINE wad_lmt  ***
1254      !!                   
1255      !! ** Purpose :   set wetting & drying mask at tracer points
1256      !!              for the current barotropic sub-step
1257      !!
1258      !! ** Method  :   ???
1259      !!
1260      !! ** Action  :  ptmsk : wetting & drying t-mask
1261      !!----------------------------------------------------------------------
1262      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
1263      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
1264      !
1265      INTEGER  ::   ji, jj   ! dummy loop indices
1266      !!----------------------------------------------------------------------
1267      !
1268      IF( ln_wd_dl_rmp ) THEN     
1269         DO_2D_11_11
1270            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
1271               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
1272               ptmsk(ji,jj) = 1._wp
1273            ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
1274               ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
1275            ELSE
1276               ptmsk(ji,jj) = 0._wp
1277            ENDIF
1278         END_2D
1279      ELSE 
1280         DO_2D_11_11
1281            IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
1282            ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
1283            ENDIF
1284         END_2D
1285      ENDIF
1286      !
1287   END SUBROUTINE wad_tmsk
1288
1289
1290   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1291      !!----------------------------------------------------------------------
1292      !!                  ***  ROUTINE wad_lmt  ***
1293      !!                   
1294      !! ** Purpose :   set wetting & drying mask at tracer points
1295      !!              for the current barotropic sub-step
1296      !!
1297      !! ** Method  :   ???
1298      !!
1299      !! ** Action  :  ptmsk : wetting & drying t-mask
1300      !!----------------------------------------------------------------------
1301      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
1302      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
1303      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
1304      !
1305      INTEGER  ::   ji, jj   ! dummy loop indices
1306      !!----------------------------------------------------------------------
1307      !
1308      DO_2D_11_10
1309         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
1310         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj) 
1311         ENDIF
1312         phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1313         pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1314      END_2D
1315      !
1316      DO_2D_10_11
1317         IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
1318         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1) 
1319         ENDIF
1320         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
1321         pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1322      END_2D
1323      !
1324   END SUBROUTINE wad_Umsk
1325
1326
1327   SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1328      !!---------------------------------------------------------------------
1329      !!                   ***  ROUTINE  wad_sp  ***
1330      !!
1331      !! ** Purpose :
1332      !!----------------------------------------------------------------------
1333      INTEGER  ::   ji ,jj               ! dummy loop indices
1334      LOGICAL  ::   ll_tmp1, ll_tmp2
1335      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
1336      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1337      !!----------------------------------------------------------------------
1338      DO_2D_00_00
1339         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
1340              &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
1341              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji+1,jj) + ht_0(ji+1,jj) )  &
1342              &                                                         > rn_wdmin1 + rn_wdmin2
1343         ll_tmp2 = ( ABS( pshn(ji+1,jj)            -  pshn(ji  ,jj))  > 1.E-12 ).AND.( &
1344              &      MAX(   pshn(ji,jj)              ,  pshn(ji+1,jj) ) >                &
1345              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1346         IF(ll_tmp1) THEN
1347            zcpx(ji,jj) = 1.0_wp
1348         ELSEIF(ll_tmp2) THEN
1349            ! no worries about  pshn(ji+1,jj) -  pshn(ji  ,jj) = 0, it won't happen ! here
1350            zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1351                 &           / (pshn(ji+1,jj) - pshn(ji  ,jj)) )
1352            zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1353         ELSE
1354            zcpx(ji,jj) = 0._wp
1355         ENDIF
1356         !
1357         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji,jj+1) ) >                &
1358              &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
1359              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji,jj+1) + ht_0(ji,jj+1) )  &
1360              &                                                       > rn_wdmin1 + rn_wdmin2
1361         ll_tmp2 = ( ABS( pshn(ji,jj)              -  pshn(ji,jj+1))  > 1.E-12 ).AND.( &
1362              &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
1363              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1364         
1365         IF(ll_tmp1) THEN
1366            zcpy(ji,jj) = 1.0_wp
1367         ELSE IF(ll_tmp2) THEN
1368            ! no worries about  pshn(ji,jj+1) -  pshn(ji,jj  ) = 0, it won't happen ! here
1369            zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1370                 &           / (pshn(ji,jj+1) - pshn(ji,jj  )) )
1371            zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
1372         ELSE
1373            zcpy(ji,jj) = 0._wp
1374         ENDIF
1375      END_2D
1376           
1377   END SUBROUTINE wad_spg
1378     
1379
1380
1381   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1382      !!----------------------------------------------------------------------
1383      !!                  ***  ROUTINE dyn_drg_init  ***
1384      !!                   
1385      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
1386      !!              the baroclinic part of the barotropic RHS
1387      !!              - compute the barotropic drag coefficients
1388      !!
1389      !! ** Method  :   computation done over the INNER domain only
1390      !!----------------------------------------------------------------------
1391      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
1392      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
1393      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
1394      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
1395      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients
1396      !
1397      INTEGER  ::   ji, jj   ! dummy loop indices
1398      INTEGER  ::   ikbu, ikbv, iktu, iktv
1399      REAL(wp) ::   zztmp
1400      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
1401      !!----------------------------------------------------------------------
1402      !
1403      !                    !==  Set the barotropic drag coef.  ==!
1404      !
1405      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
1406         
1407         DO_2D_00_00
1408            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1409            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1410         END_2D
1411      ELSE                          ! bottom friction only
1412         DO_2D_00_00
1413            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1414            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1415         END_2D
1416      ENDIF
1417      !
1418      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
1419      !
1420      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
1421         
1422         DO_2D_00_00
1423            ikbu = mbku(ji,jj)       
1424            ikbv = mbkv(ji,jj)   
1425            zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1426            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1427         END_2D
1428      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
1429         
1430         DO_2D_00_00
1431            ikbu = mbku(ji,jj)       
1432            ikbv = mbkv(ji,jj)   
1433            zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1434            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1435         END_2D
1436      ENDIF
1437      !
1438      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
1439         zztmp = -1._wp / rDt_e
1440         DO_2D_00_00
1441            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
1442                 &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
1443            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
1444                 &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
1445         END_2D
1446      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1447         
1448         DO_2D_00_00
1449            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)
1450            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)
1451         END_2D
1452      END IF
1453      !
1454      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
1455      !
1456      IF( ln_isfcav ) THEN
1457         !
1458         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
1459           
1460            DO_2D_00_00
1461               iktu = miku(ji,jj)
1462               iktv = mikv(ji,jj)
1463               zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1464               zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1465            END_2D
1466         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
1467           
1468            DO_2D_00_00
1469               iktu = miku(ji,jj)
1470               iktv = mikv(ji,jj)
1471               zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1472               zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1473            END_2D
1474         ENDIF
1475         !
1476         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1477         
1478         DO_2D_00_00
1479            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)
1480            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)
1481         END_2D
1482         !
1483      ENDIF
1484      !
1485   END SUBROUTINE dyn_drg_init
1486
1487   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
1488      &                      za0, za1, za2, za3 )   ! ==> out
1489      !!----------------------------------------------------------------------
1490      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
1491      LOGICAL ,INTENT(in   ) ::   ll_init              !
1492      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient
1493      !
1494      REAL(wp) ::   zepsilon, zgamma                   !   -      -
1495      !!----------------------------------------------------------------------
1496      !                             ! set Half-step back interpolation coefficient
1497      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
1498         za0 = 1._wp                       
1499         za1 = 0._wp                           
1500         za2 = 0._wp
1501         za3 = 0._wp
1502      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1503         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
1504         za1 =-0.1666666666666_wp                 ! za1 = gam
1505         za2 = 0.0833333333333_wp                 ! za2 = eps
1506         za3 = 0._wp             
1507      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1508         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion 
1509            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
1510            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
1511            za2 = 0.088_wp                        ! za2 = gam
1512            za3 = 0.013_wp                        ! za3 = eps
1513         ELSE                                 ! no time diffusion
1514            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1515            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1516            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1517            za1 = 1._wp - za0 - zgamma - zepsilon
1518            za2 = zgamma
1519            za3 = zepsilon
1520         ENDIF
1521      ENDIF
1522   END SUBROUTINE ts_bck_interp
1523
1524
1525   !!======================================================================
1526END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.