source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynspg_ts.F90 @ 12252

Last change on this file since 12252 was 12252, checked in by smasson, 10 months ago

rev12240_dev_r11943_MERGE_2019: same as [12251], merge trunk 12072:12248, all sette tests ok, GYRE_PISCES, AMM12, ISOMIP, VORTEX intentical to 12236

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