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_r12072_MERGE_OPTION2_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DYN/dynspg_ts.F90 @ 12251

Last change on this file since 12251 was 12251, checked in by smasson, 4 years ago

rev12232_dev_r12072_MERGE_OPTION2_2019: merge trunk 12072:12248, all sette tests ok, GYRE_PISCES, AMM12, ISOMIP, VORTEX intentical to 12210

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