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

source: NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynspg_ts.F90 @ 12077

Last change on this file since 12077 was 12077, checked in by mathiot, 3 years ago

include ENHANCE-02_ISF_nemo in UKMO merge branch

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