New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in NEMO/trunk/src/OCE/DYN – NEMO

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

Last change on this file since 15489 was 15489, checked in by jchanut, 8 months ago

Remove an useless computation in 1-way

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