source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DYN/dynspg_ts.F90 @ 11574

Last change on this file since 11574 was 11574, checked in by jchanut, 13 months ago

#2222, import changes from dev_r10973_AGRIF-01_jchanut_small_jpi_jpj (i.e. #2199)

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