New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynspg_ts.F90 @ 12107

Last change on this file since 12107 was 12107, checked in by mathiot, 3 years ago

update to the head of the trunk

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