source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DYN/dynspg_ts.F90 @ 11573

Last change on this file since 11573 was 11573, checked in by jchanut, 12 months ago

#2222, merged with trunk

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