source: NEMO/branches/2020/dev_12905_xios_restart/src/OCE/DYN/dynspg_ts.F90 @ 12950

Last change on this file since 12950 was 12950, checked in by andmirek, 6 months ago

Ticket #2462: new XIOS restart read/write interfaces

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