source: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90 @ 12326

Last change on this file since 12326 was 12326, checked in by deazer, 17 months ago

Adjusted to account for flux form advection
should now work with both types of advection and enveloping bathymetry
moved one off caluculations of depth scales to domain.F90

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