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

Last change on this file since 12095 was 12095, checked in by deazer, 10 months ago

First iteration on correcting u or v depth points when using enveloping bathy
This needs tobe done also for cases other than vec
Code is very ineffiecint as written all in dynspg_ts , this should be rewritten
to only calculate some variables at run start.

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