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

Last change on this file since 12489 was 12489, checked in by davestorkey, 8 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 78.3 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( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN    !* Read the restart file
893            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
894            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
895            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
896            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
897            IF( .NOT.ln_bt_av ) THEN
898               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
899               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
900               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
901               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
902               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
903               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
904            ENDIF
905#if defined key_agrif
906            ! Read time integrated fluxes
907            IF ( .NOT.Agrif_Root() ) THEN
908               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
909               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
910            ENDIF
911#endif
912         ELSE                                   !* Start from rest
913            IF(lwp) WRITE(numout,*)
914            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
915            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
916            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
917            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
918#if defined key_agrif
919            IF ( .NOT.Agrif_Root() ) THEN
920               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
921            ENDIF
922#endif
923         ENDIF
924         !
925      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
926         !                                   ! -------------------
927         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
928         IF( lwxios ) CALL iom_swap(      cwxios_context          )
929         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
930         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
931         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
932         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
933         !
934         IF (.NOT.ln_bt_av) THEN
935            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
936            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
937            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
938            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
939            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
940            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
941         ENDIF
942#if defined key_agrif
943         ! Save time integrated fluxes
944         IF ( .NOT.Agrif_Root() ) THEN
945            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
946            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
947         ENDIF
948#endif
949         IF( lwxios ) CALL iom_swap(      cxios_context          )
950      ENDIF
951      !
952   END SUBROUTINE ts_rst
953
954
955   SUBROUTINE dyn_spg_ts_init
956      !!---------------------------------------------------------------------
957      !!                   ***  ROUTINE dyn_spg_ts_init  ***
958      !!
959      !! ** Purpose : Set time splitting options
960      !!----------------------------------------------------------------------
961      INTEGER  ::   ji ,jj              ! dummy loop indices
962      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
963      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
964      !!----------------------------------------------------------------------
965      !
966      ! Max courant number for ext. grav. waves
967      !
968      DO_2D_11_11
969         zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
970         zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
971         zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
972      END_2D
973      !
974      zcmax = MAXVAL( zcu(:,:) )
975      CALL mpp_max( 'dynspg_ts', zcmax )
976
977      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
978      IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
979     
980      rDt_e = rn_Dt / REAL( nn_e , wp )
981      zcmax = zcmax * rDt_e
982      ! Print results
983      IF(lwp) WRITE(numout,*)
984      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
985      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
986      IF( ln_bt_auto ) THEN
987         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_e '
988         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
989      ELSE
990         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e
991      ENDIF
992
993      IF(ln_bt_av) THEN
994         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on '
995      ELSE
996         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
997      ENDIF
998      !
999      !
1000      IF(ln_bt_fw) THEN
1001         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1002      ELSE
1003         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1004      ENDIF
1005      !
1006#if defined key_agrif
1007      ! Restrict the use of Agrif to the forward case only
1008!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1009#endif
1010      !
1011      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1012      SELECT CASE ( nn_bt_flt )
1013         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1014         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
1015         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
1016         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1017      END SELECT
1018      !
1019      IF(lwp) WRITE(numout,*) ' '
1020      IF(lwp) WRITE(numout,*) '     nn_e = ', nn_e
1021      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rDt_e
1022      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1023      !
1024      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1025      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1026         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1027      ENDIF
1028      !
1029      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1030         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1031      ENDIF
1032      IF( zcmax>0.9_wp ) THEN
1033         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )         
1034      ENDIF
1035      !
1036      !                             ! Allocate time-splitting arrays
1037      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1038      !
1039      !                             ! read restart when needed
1040      CALL ts_rst( nit000, 'READ' )
1041      !
1042      IF( lwxios ) THEN
1043! define variables in restart file when writing with XIOS
1044         CALL iom_set_rstw_var_active('ub2_b')
1045         CALL iom_set_rstw_var_active('vb2_b')
1046         CALL iom_set_rstw_var_active('un_bf')
1047         CALL iom_set_rstw_var_active('vn_bf')
1048         !
1049         IF (.NOT.ln_bt_av) THEN
1050            CALL iom_set_rstw_var_active('sshbb_e')
1051            CALL iom_set_rstw_var_active('ubb_e')
1052            CALL iom_set_rstw_var_active('vbb_e')
1053            CALL iom_set_rstw_var_active('sshb_e')
1054            CALL iom_set_rstw_var_active('ub_e')
1055            CALL iom_set_rstw_var_active('vb_e')
1056         ENDIF
1057#if defined key_agrif
1058         ! Save time integrated fluxes
1059         IF ( .NOT.Agrif_Root() ) THEN
1060            CALL iom_set_rstw_var_active('ub2_i_b')
1061            CALL iom_set_rstw_var_active('vb2_i_b')
1062         ENDIF
1063#endif
1064      ENDIF
1065      !
1066   END SUBROUTINE dyn_spg_ts_init
1067
1068   
1069   SUBROUTINE dyn_cor_2D_init( Kmm )
1070      !!---------------------------------------------------------------------
1071      !!                   ***  ROUTINE dyn_cor_2D_init  ***
1072      !!
1073      !! ** Purpose : Set time splitting options
1074      !! Set arrays to remove/compute coriolis trend.
1075      !! Do it once during initialization if volume is fixed, else at each long time step.
1076      !! Note that these arrays are also used during barotropic loop. These are however frozen
1077      !! although they should be updated in the variable volume case. Not a big approximation.
1078      !! To remove this approximation, copy lines below inside barotropic loop
1079      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1080      !!
1081      !! Compute zwz = f / ( height of the water colomn )
1082      !!----------------------------------------------------------------------
1083      INTEGER,  INTENT(in)         ::  Kmm  ! Time index
1084      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
1085      REAL(wp) ::   z1_ht
1086      REAL(wp), DIMENSION(jpi,jpj) :: zhf
1087      !!----------------------------------------------------------------------
1088      !
1089      SELECT CASE( nvor_scheme )
1090      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
1091         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
1092         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
1093            DO_2D_10_10
1094               zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
1095                    &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp 
1096               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1097            END_2D
1098         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
1099            DO_2D_10_10
1100               zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
1101                    &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
1102                    &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
1103                    &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
1104               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1105            END_2D
1106         END SELECT
1107         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1108         !
1109         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1110         DO_2D_01_01
1111            ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1112            ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1113            ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1114            ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1115         END_2D
1116         !
1117      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
1118         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1119         DO_2D_01_01
1120            z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1121            ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
1122            ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
1123            ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1124            ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
1125         END_2D
1126         !
1127      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
1128         !
1129         zwz(:,:) = 0._wp
1130         zhf(:,:) = 0._wp
1131         
1132         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1133!!gm    A priori a better value should be something like :
1134!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1135!!gm                     divided by the sum of the corresponding mask
1136!!gm
1137!!           
1138         IF( .NOT.ln_sco ) THEN
1139 
1140   !!gm  agree the JC comment  : this should be done in a much clear way
1141 
1142   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1143   !     Set it to zero for the time being
1144   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
1145   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
1146   !              ENDIF
1147   !              zhf(:,:) = gdepw_0(:,:,jk+1)
1148            !
1149         ELSE
1150            !
1151            !zhf(:,:) = hbatf(:,:)
1152            DO_2D_10_10
1153               zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
1154                    &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
1155                    &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
1156                    &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
1157            END_2D
1158         ENDIF
1159         !
1160         DO jj = 1, jpjm1
1161            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1162         END DO
1163         !
1164         DO jk = 1, jpkm1
1165            DO jj = 1, jpjm1
1166               zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1167            END DO
1168         END DO
1169         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1170         ! JC: TBC. hf should be greater than 0
1171         DO_2D_11_11
1172            IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
1173         END_2D
1174         zwz(:,:) = ff_f(:,:) * zwz(:,:)
1175      END SELECT
1176     
1177   END SUBROUTINE dyn_cor_2d_init
1178
1179
1180
1181   SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
1182      !!---------------------------------------------------------------------
1183      !!                   ***  ROUTINE dyn_cor_2d  ***
1184      !!
1185      !! ** Purpose : Compute u and v coriolis trends
1186      !!----------------------------------------------------------------------
1187      INTEGER  ::   ji ,jj                             ! dummy loop indices
1188      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
1189      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phu, phv, punb, pvnb, zhU, zhV
1190      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
1191      !!----------------------------------------------------------------------
1192      SELECT CASE( nvor_scheme )
1193      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
1194         DO_2D_00_00
1195            z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1196            z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1197            zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
1198               &               * (  e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
1199               &                  + e1e2t(ji  ,jj)*ht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
1200               !
1201            zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1202               &               * (  e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
1203               &                  + e1e2t(ji,jj  )*ht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
1204         END_2D
1205         !         
1206      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
1207         DO_2D_00_00
1208            zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1209            zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1210            zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1211            zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1212            ! energy conserving formulation for planetary vorticity term
1213            zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1214            zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1215         END_2D
1216         !
1217      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
1218         DO_2D_00_00
1219            zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
1220              &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1221            zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
1222              &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1223            zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1224            zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1225         END_2D
1226         !
1227      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
1228         DO_2D_00_00
1229            zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
1230             &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
1231             &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
1232             &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1233            zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1234             &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
1235             &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
1236             &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
1237         END_2D
1238         !
1239      END SELECT
1240      !
1241   END SUBROUTINE dyn_cor_2D
1242
1243
1244   SUBROUTINE wad_tmsk( pssh, ptmsk )
1245      !!----------------------------------------------------------------------
1246      !!                  ***  ROUTINE wad_lmt  ***
1247      !!                   
1248      !! ** Purpose :   set wetting & drying mask at tracer points
1249      !!              for the current barotropic sub-step
1250      !!
1251      !! ** Method  :   ???
1252      !!
1253      !! ** Action  :  ptmsk : wetting & drying t-mask
1254      !!----------------------------------------------------------------------
1255      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
1256      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
1257      !
1258      INTEGER  ::   ji, jj   ! dummy loop indices
1259      !!----------------------------------------------------------------------
1260      !
1261      IF( ln_wd_dl_rmp ) THEN     
1262         DO_2D_11_11
1263            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
1264               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
1265               ptmsk(ji,jj) = 1._wp
1266            ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
1267               ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
1268            ELSE
1269               ptmsk(ji,jj) = 0._wp
1270            ENDIF
1271         END_2D
1272      ELSE 
1273         DO_2D_11_11
1274            IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
1275            ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
1276            ENDIF
1277         END_2D
1278      ENDIF
1279      !
1280   END SUBROUTINE wad_tmsk
1281
1282
1283   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1284      !!----------------------------------------------------------------------
1285      !!                  ***  ROUTINE wad_lmt  ***
1286      !!                   
1287      !! ** Purpose :   set wetting & drying mask at tracer points
1288      !!              for the current barotropic sub-step
1289      !!
1290      !! ** Method  :   ???
1291      !!
1292      !! ** Action  :  ptmsk : wetting & drying t-mask
1293      !!----------------------------------------------------------------------
1294      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
1295      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
1296      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
1297      !
1298      INTEGER  ::   ji, jj   ! dummy loop indices
1299      !!----------------------------------------------------------------------
1300      !
1301      DO_2D_11_10
1302         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
1303         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj) 
1304         ENDIF
1305         phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1306         pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1307      END_2D
1308      !
1309      DO_2D_10_11
1310         IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
1311         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1) 
1312         ENDIF
1313         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
1314         pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1315      END_2D
1316      !
1317   END SUBROUTINE wad_Umsk
1318
1319
1320   SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1321      !!---------------------------------------------------------------------
1322      !!                   ***  ROUTINE  wad_sp  ***
1323      !!
1324      !! ** Purpose :
1325      !!----------------------------------------------------------------------
1326      INTEGER  ::   ji ,jj               ! dummy loop indices
1327      LOGICAL  ::   ll_tmp1, ll_tmp2
1328      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
1329      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1330      !!----------------------------------------------------------------------
1331      DO_2D_00_00
1332         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
1333              &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
1334              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji+1,jj) + ht_0(ji+1,jj) )  &
1335              &                                                         > rn_wdmin1 + rn_wdmin2
1336         ll_tmp2 = ( ABS( pshn(ji+1,jj)            -  pshn(ji  ,jj))  > 1.E-12 ).AND.( &
1337              &      MAX(   pshn(ji,jj)              ,  pshn(ji+1,jj) ) >                &
1338              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1339         IF(ll_tmp1) THEN
1340            zcpx(ji,jj) = 1.0_wp
1341         ELSEIF(ll_tmp2) THEN
1342            ! no worries about  pshn(ji+1,jj) -  pshn(ji  ,jj) = 0, it won't happen ! here
1343            zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1344                 &           / (pshn(ji+1,jj) - pshn(ji  ,jj)) )
1345            zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1346         ELSE
1347            zcpx(ji,jj) = 0._wp
1348         ENDIF
1349         !
1350         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji,jj+1) ) >                &
1351              &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
1352              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji,jj+1) + ht_0(ji,jj+1) )  &
1353              &                                                       > rn_wdmin1 + rn_wdmin2
1354         ll_tmp2 = ( ABS( pshn(ji,jj)              -  pshn(ji,jj+1))  > 1.E-12 ).AND.( &
1355              &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
1356              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1357         
1358         IF(ll_tmp1) THEN
1359            zcpy(ji,jj) = 1.0_wp
1360         ELSE IF(ll_tmp2) THEN
1361            ! no worries about  pshn(ji,jj+1) -  pshn(ji,jj  ) = 0, it won't happen ! here
1362            zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1363                 &           / (pshn(ji,jj+1) - pshn(ji,jj  )) )
1364            zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
1365         ELSE
1366            zcpy(ji,jj) = 0._wp
1367         ENDIF
1368      END_2D
1369           
1370   END SUBROUTINE wad_spg
1371     
1372
1373
1374   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1375      !!----------------------------------------------------------------------
1376      !!                  ***  ROUTINE dyn_drg_init  ***
1377      !!                   
1378      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
1379      !!              the baroclinic part of the barotropic RHS
1380      !!              - compute the barotropic drag coefficients
1381      !!
1382      !! ** Method  :   computation done over the INNER domain only
1383      !!----------------------------------------------------------------------
1384      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
1385      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
1386      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
1387      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
1388      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients
1389      !
1390      INTEGER  ::   ji, jj   ! dummy loop indices
1391      INTEGER  ::   ikbu, ikbv, iktu, iktv
1392      REAL(wp) ::   zztmp
1393      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
1394      !!----------------------------------------------------------------------
1395      !
1396      !                    !==  Set the barotropic drag coef.  ==!
1397      !
1398      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
1399         
1400         DO_2D_00_00
1401            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1402            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1403         END_2D
1404      ELSE                          ! bottom friction only
1405         DO_2D_00_00
1406            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1407            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1408         END_2D
1409      ENDIF
1410      !
1411      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
1412      !
1413      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
1414         
1415         DO_2D_00_00
1416            ikbu = mbku(ji,jj)       
1417            ikbv = mbkv(ji,jj)   
1418            zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1419            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1420         END_2D
1421      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
1422         
1423         DO_2D_00_00
1424            ikbu = mbku(ji,jj)       
1425            ikbv = mbkv(ji,jj)   
1426            zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1427            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1428         END_2D
1429      ENDIF
1430      !
1431      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
1432         zztmp = -1._wp / rDt_e
1433         DO_2D_00_00
1434            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
1435                 &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
1436            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
1437                 &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
1438         END_2D
1439      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1440         
1441         DO_2D_00_00
1442            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)
1443            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)
1444         END_2D
1445      END IF
1446      !
1447      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
1448      !
1449      IF( ln_isfcav ) THEN
1450         !
1451         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
1452           
1453            DO_2D_00_00
1454               iktu = miku(ji,jj)
1455               iktv = mikv(ji,jj)
1456               zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1457               zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1458            END_2D
1459         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
1460           
1461            DO_2D_00_00
1462               iktu = miku(ji,jj)
1463               iktv = mikv(ji,jj)
1464               zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1465               zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1466            END_2D
1467         ENDIF
1468         !
1469         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1470         
1471         DO_2D_00_00
1472            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)
1473            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)
1474         END_2D
1475         !
1476      ENDIF
1477      !
1478   END SUBROUTINE dyn_drg_init
1479
1480   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
1481      &                      za0, za1, za2, za3 )   ! ==> out
1482      !!----------------------------------------------------------------------
1483      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
1484      LOGICAL ,INTENT(in   ) ::   ll_init              !
1485      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient
1486      !
1487      REAL(wp) ::   zepsilon, zgamma                   !   -      -
1488      !!----------------------------------------------------------------------
1489      !                             ! set Half-step back interpolation coefficient
1490      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
1491         za0 = 1._wp                       
1492         za1 = 0._wp                           
1493         za2 = 0._wp
1494         za3 = 0._wp
1495      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1496         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
1497         za1 =-0.1666666666666_wp                 ! za1 = gam
1498         za2 = 0.0833333333333_wp                 ! za2 = eps
1499         za3 = 0._wp             
1500      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1501         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion 
1502            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
1503            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
1504            za2 = 0.088_wp                        ! za2 = gam
1505            za3 = 0.013_wp                        ! za3 = eps
1506         ELSE                                 ! no time diffusion
1507            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1508            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1509            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1510            za1 = 1._wp - za0 - zgamma - zepsilon
1511            za2 = zgamma
1512            za3 = zepsilon
1513         ENDIF
1514      ENDIF
1515   END SUBROUTINE ts_bck_interp
1516
1517
1518   !!======================================================================
1519END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.