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/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90 @ 11234

Last change on this file since 11234 was 11234, checked in by girrmann, 5 years ago

dev_r10984_HPC-13 : cleaning of dyn_spg_ts routine, does not change the results, see #2285

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