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

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

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

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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