source: NEMO/branches/2019/dev_r11756_SI3restart_XIOS/src/OCE/DYN/dynspg_ts.F90 @ 11837

Last change on this file since 11837 was 11837, checked in by andmirek, 18 months ago

ticket #2323 read SI3 restart with XIOS

  • Property svn:keywords set to Id
File size: 81.5 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(lrxios) CALL iom_swap(crxios_context)
937         IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN    !* Read the restart file
938            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
939            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
940            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
941            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
942            IF( .NOT.ln_bt_av ) THEN
943               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
944               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
945               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
946               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
947               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
948               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
949            ENDIF
950#if defined key_agrif
951            ! Read time integrated fluxes
952            IF ( .NOT.Agrif_Root() ) THEN
953               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
954               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
955            ENDIF
956#endif
957         IF(lrxios) CALL iom_swap(cxios_context)
958         !
959         ELSE                                   !* Start from rest
960            IF(lwp) WRITE(numout,*)
961            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
962            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
963            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
964            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
965#if defined key_agrif
966            IF ( .NOT.Agrif_Root() ) THEN
967               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
968            ENDIF
969#endif
970         ENDIF
971         !
972      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
973         !                                   ! -------------------
974         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
975         IF( lwxios ) CALL iom_swap(      cwxios_context          )
976         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
977         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
978         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
979         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
980         !
981         IF (.NOT.ln_bt_av) THEN
982            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
983            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
984            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
985            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
986            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
987            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
988         ENDIF
989#if defined key_agrif
990         ! Save time integrated fluxes
991         IF ( .NOT.Agrif_Root() ) THEN
992            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
993            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
994         ENDIF
995#endif
996         IF( lwxios ) CALL iom_swap(      cxios_context          )
997      ENDIF
998      !
999   END SUBROUTINE ts_rst
1000
1001
1002   SUBROUTINE dyn_spg_ts_init
1003      !!---------------------------------------------------------------------
1004      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1005      !!
1006      !! ** Purpose : Set time splitting options
1007      !!----------------------------------------------------------------------
1008      INTEGER  ::   ji ,jj              ! dummy loop indices
1009      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1010      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
1011      INTEGER  :: inum
1012      !!----------------------------------------------------------------------
1013      !
1014      ! Max courant number for ext. grav. waves
1015      !
1016      DO jj = 1, jpj
1017         DO ji =1, jpi
1018            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1019            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1020            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1021         END DO
1022      END DO
1023      !
1024      zcmax = MAXVAL( zcu(:,:) )
1025      CALL mpp_max( 'dynspg_ts', zcmax )
1026
1027      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1028      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1029     
1030      rdtbt = rdt / REAL( nn_baro , wp )
1031      zcmax = zcmax * rdtbt
1032      ! Print results
1033      IF(lwp) WRITE(numout,*)
1034      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1035      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1036      IF( ln_bt_auto ) THEN
1037         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
1038         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1039      ELSE
1040         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
1041      ENDIF
1042
1043      IF(ln_bt_av) THEN
1044         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
1045      ELSE
1046         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1047      ENDIF
1048      !
1049      !
1050      IF(ln_bt_fw) THEN
1051         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1052      ELSE
1053         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1054      ENDIF
1055      !
1056#if defined key_agrif
1057      ! Restrict the use of Agrif to the forward case only
1058!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1059#endif
1060      !
1061      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1062      SELECT CASE ( nn_bt_flt )
1063         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1064         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1065         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1066         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1067      END SELECT
1068      !
1069      IF(lwp) WRITE(numout,*) ' '
1070      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1071      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1072      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1073      !
1074      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1075      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1076         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1077      ENDIF
1078      !
1079      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1080         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1081      ENDIF
1082      IF( zcmax>0.9_wp ) THEN
1083         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1084      ENDIF
1085      !
1086      !                             ! Allocate time-splitting arrays
1087      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1088      !
1089      !                             ! read restart when needed
1090      CALL ts_rst( nit000, 'READ' )
1091      !
1092      IF( lwxios ) THEN
1093! define variables in restart file when writing with XIOS
1094         CALL iom_set_rstw_var_active(rst_wfields, 'ub2_b')
1095         CALL iom_set_rstw_var_active(rst_wfields, 'vb2_b')
1096         CALL iom_set_rstw_var_active(rst_wfields, 'un_bf')
1097         CALL iom_set_rstw_var_active(rst_wfields, 'vn_bf')
1098         !
1099         IF (.NOT.ln_bt_av) THEN
1100            CALL iom_set_rstw_var_active(rst_wfields, 'sshbb_e')
1101            CALL iom_set_rstw_var_active(rst_wfields, 'ubb_e')
1102            CALL iom_set_rstw_var_active(rst_wfields, 'vbb_e')
1103            CALL iom_set_rstw_var_active(rst_wfields, 'sshb_e')
1104            CALL iom_set_rstw_var_active(rst_wfields, 'ub_e')
1105            CALL iom_set_rstw_var_active(rst_wfields, 'vb_e')
1106         ENDIF
1107#if defined key_agrif
1108         ! Save time integrated fluxes
1109         IF ( .NOT.Agrif_Root() ) THEN
1110            CALL iom_set_rstw_var_active(rst_wfields, 'ub2_i_b')
1111            CALL iom_set_rstw_var_active(rst_wfields, 'vb2_i_b')
1112         ENDIF
1113#endif
1114      ENDIF
1115      !
1116   END SUBROUTINE dyn_spg_ts_init
1117
1118   
1119   SUBROUTINE dyn_cor_2d_init
1120      !!---------------------------------------------------------------------
1121      !!                   ***  ROUTINE dyn_cor_2d_init  ***
1122      !!
1123      !! ** Purpose : Set time splitting options
1124      !! Set arrays to remove/compute coriolis trend.
1125      !! Do it once during initialization if volume is fixed, else at each long time step.
1126      !! Note that these arrays are also used during barotropic loop. These are however frozen
1127      !! although they should be updated in the variable volume case. Not a big approximation.
1128      !! To remove this approximation, copy lines below inside barotropic loop
1129      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1130      !!
1131      !! Compute zwz = f / ( height of the water colomn )
1132      !!----------------------------------------------------------------------
1133      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
1134      REAL(wp) ::   z1_ht
1135      REAL(wp), DIMENSION(jpi,jpj) :: zhf
1136      !!----------------------------------------------------------------------
1137      !
1138      SELECT CASE( nvor_scheme )
1139      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
1140         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
1141         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
1142            DO jj = 1, jpjm1
1143               DO ji = 1, jpim1
1144                  zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
1145                       &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp 
1146                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1147               END DO
1148            END DO
1149         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
1150            DO jj = 1, jpjm1
1151               DO ji = 1, jpim1
1152                  zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
1153                       &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
1154                       &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
1155                       &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
1156                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1157               END DO
1158            END DO
1159         END SELECT
1160         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1161         !
1162         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1163         DO jj = 2, jpj
1164            DO ji = 2, jpi
1165               ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1166               ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1167               ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1168               ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1169            END DO
1170         END DO
1171         !
1172      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
1173         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1174         DO jj = 2, jpj
1175            DO ji = 2, jpi
1176               z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
1177               ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
1178               ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
1179               ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1180               ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
1181            END DO
1182         END DO
1183         !
1184      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
1185         !
1186         zwz(:,:) = 0._wp
1187         zhf(:,:) = 0._wp
1188         
1189         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1190!!gm    A priori a better value should be something like :
1191!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1192!!gm                     divided by the sum of the corresponding mask
1193!!gm
1194!!           
1195         IF( .NOT.ln_sco ) THEN
1196 
1197   !!gm  agree the JC comment  : this should be done in a much clear way
1198 
1199   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1200   !     Set it to zero for the time being
1201   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
1202   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
1203   !              ENDIF
1204   !              zhf(:,:) = gdepw_0(:,:,jk+1)
1205            !
1206         ELSE
1207            !
1208            !zhf(:,:) = hbatf(:,:)
1209            DO jj = 1, jpjm1
1210               DO ji = 1, jpim1
1211                  zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
1212                       &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
1213                       &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
1214                       &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
1215               END DO
1216            END DO
1217         ENDIF
1218         !
1219         DO jj = 1, jpjm1
1220            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1221         END DO
1222         !
1223         DO jk = 1, jpkm1
1224            DO jj = 1, jpjm1
1225               zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1226            END DO
1227         END DO
1228         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1229         ! JC: TBC. hf should be greater than 0
1230         DO jj = 1, jpj
1231            DO ji = 1, jpi
1232               IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
1233            END DO
1234         END DO
1235         zwz(:,:) = ff_f(:,:) * zwz(:,:)
1236      END SELECT
1237     
1238   END SUBROUTINE dyn_cor_2d_init
1239
1240
1241
1242   SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV,    zu_trd, zv_trd   )
1243      !!---------------------------------------------------------------------
1244      !!                   ***  ROUTINE dyn_cor_2d  ***
1245      !!
1246      !! ** Purpose : Compute u and v coriolis trends
1247      !!----------------------------------------------------------------------
1248      INTEGER  ::   ji ,jj                             ! dummy loop indices
1249      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
1250      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV
1251      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
1252      !!----------------------------------------------------------------------
1253      SELECT CASE( nvor_scheme )
1254      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
1255         DO jj = 2, jpjm1
1256            DO ji = 2, jpim1
1257               z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )
1258               z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )
1259               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
1260                  &               * (  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) )   &
1261                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
1262                  !
1263               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1264                  &               * (  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) )   & 
1265                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
1266            END DO 
1267         END DO 
1268         !         
1269      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
1270         DO jj = 2, jpjm1
1271            DO ji = fs_2, fs_jpim1   ! vector opt.
1272               zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1273               zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1274               zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1275               zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1276               ! energy conserving formulation for planetary vorticity term
1277               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1278               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1279            END DO
1280         END DO
1281         !
1282      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
1283         DO jj = 2, jpjm1
1284            DO ji = fs_2, fs_jpim1   ! vector opt.
1285               zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
1286                 &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1287               zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
1288                 &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1289               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1290               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1291            END DO
1292         END DO
1293         !
1294      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
1295         DO jj = 2, jpjm1
1296            DO ji = fs_2, fs_jpim1   ! vector opt.
1297               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
1298                &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
1299                &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
1300                &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1301               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1302                &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
1303                &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
1304                &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
1305            END DO
1306         END DO
1307         !
1308      END SELECT
1309      !
1310   END SUBROUTINE dyn_cor_2D
1311
1312
1313   SUBROUTINE wad_tmsk( pssh, ptmsk )
1314      !!----------------------------------------------------------------------
1315      !!                  ***  ROUTINE wad_lmt  ***
1316      !!                   
1317      !! ** Purpose :   set wetting & drying mask at tracer points
1318      !!              for the current barotropic sub-step
1319      !!
1320      !! ** Method  :   ???
1321      !!
1322      !! ** Action  :  ptmsk : wetting & drying t-mask
1323      !!----------------------------------------------------------------------
1324      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
1325      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
1326      !
1327      INTEGER  ::   ji, jj   ! dummy loop indices
1328      !!----------------------------------------------------------------------
1329      !
1330      IF( ln_wd_dl_rmp ) THEN     
1331         DO jj = 1, jpj
1332            DO ji = 1, jpi                   
1333               IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
1334                  !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
1335                  ptmsk(ji,jj) = 1._wp
1336               ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
1337                  ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
1338               ELSE
1339                  ptmsk(ji,jj) = 0._wp
1340               ENDIF
1341            END DO
1342         END DO
1343      ELSE 
1344         DO jj = 1, jpj
1345            DO ji = 1, jpi                             
1346               IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
1347               ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
1348               ENDIF
1349            END DO
1350         END DO
1351      ENDIF
1352      !
1353   END SUBROUTINE wad_tmsk
1354
1355
1356   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1357      !!----------------------------------------------------------------------
1358      !!                  ***  ROUTINE wad_lmt  ***
1359      !!                   
1360      !! ** Purpose :   set wetting & drying mask at tracer points
1361      !!              for the current barotropic sub-step
1362      !!
1363      !! ** Method  :   ???
1364      !!
1365      !! ** Action  :  ptmsk : wetting & drying t-mask
1366      !!----------------------------------------------------------------------
1367      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
1368      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
1369      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
1370      !
1371      INTEGER  ::   ji, jj   ! dummy loop indices
1372      !!----------------------------------------------------------------------
1373      !
1374      DO jj = 1, jpj
1375         DO ji = 1, jpim1   ! not jpi-column
1376            IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
1377            ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj) 
1378            ENDIF
1379            phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1380            pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1381         END DO
1382      END DO
1383      !
1384      DO jj = 1, jpjm1   ! not jpj-row
1385         DO ji = 1, jpi
1386            IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
1387            ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1) 
1388            ENDIF
1389            phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
1390            pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1391         END DO
1392      END DO
1393      !
1394   END SUBROUTINE wad_Umsk
1395
1396
1397   SUBROUTINE wad_spg( sshn, zcpx, zcpy )
1398      !!---------------------------------------------------------------------
1399      !!                   ***  ROUTINE  wad_sp  ***
1400      !!
1401      !! ** Purpose :
1402      !!----------------------------------------------------------------------
1403      INTEGER  ::   ji ,jj               ! dummy loop indices
1404      LOGICAL  ::   ll_tmp1, ll_tmp2
1405      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: sshn
1406      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1407      !!----------------------------------------------------------------------
1408      DO jj = 2, jpjm1
1409         DO ji = 2, jpim1 
1410            ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
1411                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
1412                 &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
1413                 &                                                         > rn_wdmin1 + rn_wdmin2
1414            ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
1415                 &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
1416                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1417            IF(ll_tmp1) THEN
1418               zcpx(ji,jj) = 1.0_wp
1419            ELSEIF(ll_tmp2) THEN
1420               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
1421               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
1422                    &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
1423               zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1424            ELSE
1425               zcpx(ji,jj) = 0._wp
1426            ENDIF
1427            !
1428            ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
1429                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
1430                 &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
1431                 &                                                       > rn_wdmin1 + rn_wdmin2
1432            ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
1433                 &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
1434                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1435           
1436            IF(ll_tmp1) THEN
1437               zcpy(ji,jj) = 1.0_wp
1438            ELSE IF(ll_tmp2) THEN
1439               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
1440               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
1441                    &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
1442               zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
1443            ELSE
1444               zcpy(ji,jj) = 0._wp
1445            ENDIF
1446         END DO
1447      END DO
1448           
1449   END SUBROUTINE wad_spg
1450     
1451
1452
1453   SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1454      !!----------------------------------------------------------------------
1455      !!                  ***  ROUTINE dyn_drg_init  ***
1456      !!                   
1457      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
1458      !!              the baroclinic part of the barotropic RHS
1459      !!              - compute the barotropic drag coefficients
1460      !!
1461      !! ** Method  :   computation done over the INNER domain only
1462      !!----------------------------------------------------------------------
1463      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
1464      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pCdU_u , pCdU_v    ! barotropic drag coefficients
1465      !
1466      INTEGER  ::   ji, jj   ! dummy loop indices
1467      INTEGER  ::   ikbu, ikbv, iktu, iktv
1468      REAL(wp) ::   zztmp
1469      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
1470      !!----------------------------------------------------------------------
1471      !
1472      !                    !==  Set the barotropic drag coef.  ==!
1473      !
1474      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
1475         
1476         DO jj = 2, jpjm1
1477            DO ji = 2, jpim1     ! INNER domain
1478               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1479               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1480            END DO
1481         END DO
1482      ELSE                          ! bottom friction only
1483         DO jj = 2, jpjm1
1484            DO ji = 2, jpim1  ! INNER domain
1485               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1486               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1487            END DO
1488         END DO
1489      ENDIF
1490      !
1491      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
1492      !
1493      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
1494         
1495         DO jj = 2, jpjm1
1496            DO ji = 2, jpim1  ! INNER domain
1497               ikbu = mbku(ji,jj)       
1498               ikbv = mbkv(ji,jj)   
1499               zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj)
1500               zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
1501            END DO
1502         END DO
1503      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
1504         
1505         DO jj = 2, jpjm1
1506            DO ji = 2, jpim1   ! INNER domain
1507               ikbu = mbku(ji,jj)       
1508               ikbv = mbkv(ji,jj)   
1509               zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj)
1510               zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
1511            END DO
1512         END DO
1513      ENDIF
1514      !
1515      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
1516         zztmp = -1._wp / rdtbt
1517         DO jj = 2, jpjm1
1518            DO ji = 2, jpim1    ! INNER domain
1519               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
1520                    &                              r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
1521               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
1522                    &                              r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
1523            END DO
1524         END DO
1525      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1526         
1527         DO jj = 2, jpjm1
1528            DO ji = 2, jpim1    ! INNER domain
1529               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)
1530               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)
1531            END DO
1532         END DO
1533      END IF
1534      !
1535      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
1536      !
1537      IF( ln_isfcav ) THEN
1538         !
1539         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
1540           
1541            DO jj = 2, jpjm1
1542               DO ji = 2, jpim1   ! INNER domain
1543                  iktu = miku(ji,jj)
1544                  iktv = mikv(ji,jj)
1545                  zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj)
1546                  zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
1547               END DO
1548            END DO
1549         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
1550           
1551            DO jj = 2, jpjm1
1552               DO ji = 2, jpim1      ! INNER domain
1553                  iktu = miku(ji,jj)
1554                  iktv = mikv(ji,jj)
1555                  zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj)
1556                  zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
1557               END DO
1558            END DO
1559         ENDIF
1560         !
1561         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1562         
1563         DO jj = 2, jpjm1
1564            DO ji = 2, jpim1    ! INNER domain
1565               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)
1566               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)
1567            END DO
1568         END DO
1569         !
1570      ENDIF
1571      !
1572   END SUBROUTINE dyn_drg_init
1573
1574   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
1575      &                      za0, za1, za2, za3 )   ! ==> out
1576      !!----------------------------------------------------------------------
1577      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
1578      LOGICAL ,INTENT(in   ) ::   ll_init              !
1579      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient
1580      !
1581      REAL(wp) ::   zepsilon, zgamma                   !   -      -
1582      !!----------------------------------------------------------------------
1583      !                             ! set Half-step back interpolation coefficient
1584      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
1585         za0 = 1._wp                       
1586         za1 = 0._wp                           
1587         za2 = 0._wp
1588         za3 = 0._wp
1589      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1590         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
1591         za1 =-0.1666666666666_wp                 ! za1 = gam
1592         za2 = 0.0833333333333_wp                 ! za2 = eps
1593         za3 = 0._wp             
1594      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1595         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion 
1596            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
1597            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
1598            za2 = 0.088_wp                        ! za2 = gam
1599            za3 = 0.013_wp                        ! za3 = eps
1600         ELSE                                 ! no time diffusion
1601            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1602            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1603            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1604            za1 = 1._wp - za0 - zgamma - zepsilon
1605            za2 = zgamma
1606            za3 = zepsilon
1607         ENDIF
1608      ENDIF
1609   END SUBROUTINE ts_bck_interp
1610
1611
1612   !!======================================================================
1613END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.