source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynspg_ts.F90 @ 12250

Last change on this file since 12250 was 12250, checked in by smasson, 11 months ago

rev12240_dev_r11943_MERGE_2019: same as [12246], add modifications from dev_r12114_ticket_2263, results unchanged except SPITZ12 as explained in #2263

  • Property svn:keywords set to Id
File size: 82.4 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 tide_mod        !
47   USE sbcwave         ! surface wave
48   USE diatmb          ! Top,middle,bottom output
49#if defined key_agrif
50   USE agrif_oce_interp ! agrif
51   USE agrif_oce
52#endif
53#if defined key_asminc   
54   USE asminc          ! Assimilation increment
55#endif
56   !
57   USE in_out_manager  ! I/O manager
58   USE lib_mpp         ! distributed memory computing library
59   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl          ! Print control
61   USE iom             ! IOM library
62   USE restart         ! only for lrst_oce
63   USE diatmb          ! Top,middle,bottom output
64
65   USE iom   ! to remove
66
67   IMPLICIT NONE
68   PRIVATE
69
70   PUBLIC dyn_spg_ts        ! called by dyn_spg
71   PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
72
73   !! Time filtered arrays at baroclinic time step:
74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step
75   !
76   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
77   REAL(wp),SAVE :: rdtbt       ! Barotropic time step
78   !
79   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields
80   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwz                ! ff_f/h at F points
81   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftnw, ftne         ! triad of coriolis parameter
82   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftsw, ftse         ! (only used with een vorticity scheme)
83
84   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios
85   REAL(wp) ::   r1_8  = 0.125_wp         !
86   REAL(wp) ::   r1_4  = 0.25_wp          !
87   REAL(wp) ::   r1_2  = 0.5_wp           !
88
89   !! * Substitutions
90#  include "vectopt_loop_substitute.h90"
91   !!----------------------------------------------------------------------
92   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
93   !! $Id$
94   !! Software governed by the CeCILL license (see ./LICENSE)
95   !!----------------------------------------------------------------------
96CONTAINS
97
98   INTEGER FUNCTION dyn_spg_ts_alloc()
99      !!----------------------------------------------------------------------
100      !!                  ***  routine dyn_spg_ts_alloc  ***
101      !!----------------------------------------------------------------------
102      INTEGER :: ierr(3)
103      !!----------------------------------------------------------------------
104      ierr(:) = 0
105      !
106      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
107      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
108         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2)   )
109         !
110      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
111      !
112      dyn_spg_ts_alloc = MAXVAL( ierr(:) )
113      !
114      CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
115      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
116      !
117   END FUNCTION dyn_spg_ts_alloc
118
119
120   SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
121      !!----------------------------------------------------------------------
122      !!
123      !! ** Purpose : - Compute the now trend due to the explicit time stepping
124      !!              of the quasi-linear barotropic system, and add it to the
125      !!              general momentum trend.
126      !!
127      !! ** Method  : - split-explicit schem (time splitting) :
128      !!      Barotropic variables are advanced from internal time steps
129      !!      "n"   to "n+1" if ln_bt_fw=T
130      !!      or from
131      !!      "n-1" to "n+1" if ln_bt_fw=F
132      !!      thanks to a generalized forward-backward time stepping (see ref. below).
133      !!
134      !! ** Action :
135      !!      -Update the filtered free surface at step "n+1"      : pssh(:,:,Kaa)
136      !!      -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa)
137      !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
138      !!      These are used to advect tracers and are compliant with discrete
139      !!      continuity equation taken at the baroclinic time steps. This
140      !!      ensures tracers conservation.
141      !!      - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component.
142      !!
143      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
144      !!---------------------------------------------------------------------
145      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
146      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
147      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
148      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
149      !
150      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
151      LOGICAL  ::   ll_fw_start           ! =T : forward integration
152      LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
153      INTEGER  ::   noffset               ! local integers  : time offset for bdy update
154      REAL(wp) ::   r1_2dt_b, z1_hu, z1_hv          ! local scalars
155      REAL(wp) ::   za0, za1, za2, za3              !   -      -
156      REAL(wp) ::   zmdi, zztmp, zldg               !   -      -
157      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      -
158      REAL(wp) ::   zun_save, zvn_save              !   -      -
159      REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc
160      REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
161      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e
162      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
163      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
164      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
165      !
166      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.
167
168      INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point
169
170      REAL(wp) ::   zepsilon, zgamma            !   -      -
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy   ! Wetting/Dying gravity filter coef.
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
174      REAL(wp) ::   zt0substep !   Time of day at the beginning of the time substep
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, pt_offset= REAL(noffset+1,wp) )
446         ! Update tide potential at the beginning of current time substep
447         IF( ln_tide_pot .AND. ln_tide ) THEN
448            zt0substep = REAL(nsec_day, wp) - 0.5_wp*rdt + (jn + noffset - 1) * rdt / REAL(nn_baro, wp)
449            CALL upd_tide(zt0substep, Kmm)
450         END IF
451         !
452         !                    !==  extrapolation at mid-step  ==!   (jn+1/2)
453         !
454         !                       !* Set extrapolation coefficients for predictor step:
455         IF ((jn<3).AND.ll_init) THEN      ! Forward           
456           za1 = 1._wp                                         
457           za2 = 0._wp                       
458           za3 = 0._wp                       
459         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
460           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
461           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
462           za3 =  0.281105_wp              ! za3 = bet
463         ENDIF
464         !
465         !                       !* Extrapolate barotropic velocities at mid-step (jn+1/2)
466         !--        m+1/2               m                m-1           m-2       --!
467         !--       u      = (3/2+beta) u   -(1/2+2beta) u      + beta u          --!
468         !-------------------------------------------------------------------------!
469         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
470         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
471
472         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
473            !                                             !  ------------------
474            ! Extrapolate Sea Level at step jit+0.5:
475            !--         m+1/2                 m                  m-1             m-2       --!
476            !--      ssh      = (3/2+beta) ssh   -(1/2+2beta) ssh      + beta ssh          --!
477            !--------------------------------------------------------------------------------!
478            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
479           
480            ! set wetting & drying mask at tracer points for this barotropic mid-step
481            IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask )
482            !
483            !                          ! ocean t-depth at mid-step
484            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
485            !
486            !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
487            DO jj = 1, jpj
488               DO ji = 1, jpim1   ! not jpi-column
489                  zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
490                       &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
491                       &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
492               END DO
493            END DO
494            DO jj = 1, jpjm1        ! not jpj-row
495               DO ji = 1, jpi
496                  zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
497                       &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
498                       &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
499               END DO
500            END DO
501            !
502         ENDIF
503         !
504         !                    !==  after SSH  ==!   (jn+1)
505         !
506         !                             ! update (ua_e,va_e) to enforce volume conservation at open boundaries
507         !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
508         IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
509         !     
510         !                             ! resulting flux at mid-step (not over the full domain)
511         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
512         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
513         !
514#if defined key_agrif
515         ! Set fluxes during predictor step to ensure volume conservation
516         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
517#endif
518         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
519
520         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where
521            !                          ! the direction of the flow is from dry cells
522            CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )   ! not jpi colomn for U, not jpj row for V
523            !
524         ENDIF   
525         !
526         !
527         !     Compute Sea Level at step jit+1
528         !--           m+1        m                               m+1/2          --!
529         !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
530         !-------------------------------------------------------------------------!
531         DO jj = 2, jpjm1        ! INNER domain                             
532            DO ji = 2, jpim1
533               zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
534               ssha_e(ji,jj) = (  sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
535            END DO
536         END DO
537         !
538         CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp )
539         !
540         !                             ! Sum over sub-time-steps to compute advective velocities
541         za2 = wgtbtp2(jn)             ! zhU, zhV hold fluxes extrapolated at jn+0.5
542         un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
543         vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
544         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
545         IF ( ln_wd_dl_bc ) THEN
546            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
547            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row
548         END IF
549         !
550         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
551         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
552#if defined key_agrif
553         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
554#endif
555         
556         ! Sea Surface Height at u-,v-points (vvl case only)
557         IF( .NOT.ln_linssh ) THEN                               
558            DO jj = 2, jpjm1   ! INNER domain, will be extended to whole domain later
559               DO ji = 2, jpim1      ! NO Vector Opt.
560                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
561                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
562                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
563                  zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
564                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
565                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
566               END DO
567            END DO
568         ENDIF   
569         !         
570         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
571         !--            m+1/2           m+1              m               m-1              m-2     --!
572         !--        ssh'    =  za0 * ssh     +  za1 * ssh   +  za2 * ssh      +  za3 * ssh        --!
573         !------------------------------------------------------------------------------------------!
574         CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 )   ! coeficients of the interpolation
575         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
576            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
577         !
578         !                             ! Surface pressure gradient
579         zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
580         DO jj = 2, jpjm1                           
581            DO ji = 2, jpim1
582               zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
583               zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
584            END DO
585         END DO
586         IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient
587            CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters
588            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
589            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
590         ENDIF
591         !
592         ! Add Coriolis trend:
593         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
594         ! at each time step. We however keep them constant here for optimization.
595         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
596         CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
597         !
598         ! Add tidal astronomical forcing if defined
599         IF ( ln_tide .AND. ln_tide_pot ) THEN
600            DO jj = 2, jpjm1
601               DO ji = fs_2, fs_jpim1   ! vector opt.
602                  zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
603                  zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
604               END DO
605            END DO
606         ENDIF
607         !
608         ! Add bottom stresses:
609!jth do implicitly instead
610         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
611            DO jj = 2, jpjm1
612               DO ji = fs_2, fs_jpim1   ! vector opt.
613                  zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
614                  zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
615               END DO
616            END DO
617         ENDIF
618         !
619         ! Set next velocities:
620         !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn)
621         !--                              VECTOR FORM
622         !--   m+1                 m               /                                                       m+1/2           \    --!
623         !--  u     =             u   + delta_t' * \         (1-r)*g * grad_x( ssh') -         f * k vect u      +     frc /    --!
624         !--                                                                                                                    --!
625         !--                             FLUX FORM                                                                              --!
626         !--  m+1   __1__  /  m    m               /  m+1/2                             m+1/2              m+1/2    n      \ \  --!
627         !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --!
628         !--         h     \                                                                                                 /  --!
629         !------------------------------------------------------------------------------------------------------------------------!
630         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
631            DO jj = 2, jpjm1
632               DO ji = fs_2, fs_jpim1   ! vector opt.
633                  ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
634                            &     + rdtbt * (                   zu_spg(ji,jj)   &
635                            &                                 + zu_trd(ji,jj)   &
636                            &                                 + zu_frc(ji,jj) ) & 
637                            &   ) * ssumask(ji,jj)
638
639                  va_e(ji,jj) = (                                 vn_e(ji,jj)   &
640                            &     + rdtbt * (                   zv_spg(ji,jj)   &
641                            &                                 + zv_trd(ji,jj)   &
642                            &                                 + zv_frc(ji,jj) ) &
643                            &   ) * ssvmask(ji,jj)
644               END DO
645            END DO
646            !
647         ELSE                           !* Flux form
648            DO jj = 2, jpjm1
649               DO ji = 2, jpim1
650                  !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
651                  !                    ! backward interpolated depth used in spg terms at jn+1/2
652                  zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
653                       &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
654                  zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
655                       &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
656                  !                    ! inverse depth at jn+1
657                  z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
658                  z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
659                  !
660                  ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
661                       &            + rdtbt * (  zhu_bck        * zu_spg (ji,jj)  &   !
662                       &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
663                       &                       +  hu(ji,jj,Kmm) * zu_frc (ji,jj)  )   ) * z1_hu
664                  !
665                  va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      &
666                       &            + rdtbt * (  zhv_bck        * zv_spg (ji,jj)  &   !
667                       &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   !
668                       &                       +  hv(ji,jj,Kmm) * zv_frc (ji,jj)  )   ) * z1_hv
669               END DO
670            END DO
671         ENDIF
672!jth implicit bottom friction:
673         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
674            DO jj = 2, jpjm1
675               DO ji = fs_2, fs_jpim1   ! vector opt.
676                     ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
677                     va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
678               END DO
679            END DO
680         ENDIF
681       
682         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
683            hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
684            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
685            hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
686            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
687            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
688                 &                         , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
689                 &                         , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
690         ELSE
691            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
692         ENDIF
693         !
694         !
695         !                                                 ! open boundaries
696         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
697#if defined key_agrif                                                           
698         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
699#endif
700         !                                             !* Swap
701         !                                             !  ----
702         ubb_e  (:,:) = ub_e  (:,:)
703         ub_e   (:,:) = un_e  (:,:)
704         un_e   (:,:) = ua_e  (:,:)
705         !
706         vbb_e  (:,:) = vb_e  (:,:)
707         vb_e   (:,:) = vn_e  (:,:)
708         vn_e   (:,:) = va_e  (:,:)
709         !
710         sshbb_e(:,:) = sshb_e(:,:)
711         sshb_e (:,:) = sshn_e(:,:)
712         sshn_e (:,:) = ssha_e(:,:)
713
714         !                                             !* Sum over whole bt loop
715         !                                             !  ----------------------
716         za1 = wgtbtp1(jn)                                   
717         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
718            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
719            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
720         ELSE                                       ! Sum transports
721            IF ( .NOT.ln_wd_dl ) THEN 
722               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
723               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
724            ELSE
725               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
726               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
727            END IF
728         ENDIF
729         !                                          ! Sum sea level
730         pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
731
732         !                                                 ! ==================== !
733      END DO                                               !        end loop      !
734      !                                                    ! ==================== !
735      ! -----------------------------------------------------------------------------
736      ! Phase 3. update the general trend with the barotropic trend
737      ! -----------------------------------------------------------------------------
738      !
739      ! Set advection velocity correction:
740      IF (ln_bt_fw) THEN
741         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
742            DO jj = 1, jpj
743               DO ji = 1, jpi
744                  zun_save = un_adv(ji,jj)
745                  zvn_save = vn_adv(ji,jj)
746                  !                          ! apply the previously computed correction
747                  un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
748                  vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
749                  !                          ! Update corrective fluxes for next time step
750                  un_bf(ji,jj)  = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
751                  vn_bf(ji,jj)  = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
752                  !                          ! Save integrated transport for next computation
753                  ub2_b(ji,jj) = zun_save
754                  vb2_b(ji,jj) = zvn_save
755               END DO
756            END DO
757         ELSE
758            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero
759            vn_bf(:,:) = 0._wp
760            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
761            vb2_b(:,:) = vn_adv(:,:)
762         END IF
763      ENDIF
764
765
766      !
767      ! Update barotropic trend:
768      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
769         DO jk=1,jpkm1
770            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_2dt_b
771            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_2dt_b
772         END DO
773      ELSE
774         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
775         DO jj = 1, jpjm1
776            DO ji = 1, jpim1      ! NO Vector Opt.
777               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
778                  &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      &
779                  &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
780               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
781                  &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      &
782                  &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
783            END DO
784         END DO
785         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
786         !
787         DO jk=1,jpkm1
788            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_2dt_b
789            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_2dt_b
790         END DO
791         ! Save barotropic velocities not transport:
792         puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
793         pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
794      ENDIF
795
796
797      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
798      DO jk = 1, jpkm1
799         puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
800         pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)
801      END DO
802
803      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
804         DO jk = 1, jpkm1
805            puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
806                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
807            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
808                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 
809         END DO
810      END IF
811
812     
813      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
814      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
815      !
816#if defined key_agrif
817      ! Save time integrated fluxes during child grid integration
818      ! (used to update coarse grid transports at next time step)
819      !
820      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
821         IF( Agrif_NbStepint() == 0 ) THEN
822            ub2_i_b(:,:) = 0._wp
823            vb2_i_b(:,:) = 0._wp
824         END IF
825         !
826         za1 = 1._wp / REAL(Agrif_rhot(), wp)
827         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
828         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
829      ENDIF
830#endif     
831      !                                   !* write time-spliting arrays in the restart
832      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
833      !
834      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
835      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
836      !
837      IF( ln_diatmb ) THEN
838         CALL iom_put( "baro_u" , puu_b(:,:,Kmm)*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
839         CALL iom_put( "baro_v" , pvv_b(:,:,Kmm)*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
840      ENDIF
841      !
842   END SUBROUTINE dyn_spg_ts
843
844
845   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
846      !!---------------------------------------------------------------------
847      !!                   ***  ROUTINE ts_wgt  ***
848      !!
849      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
850      !!----------------------------------------------------------------------
851      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
852      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
853      INTEGER, INTENT(inout) :: jpit      ! cycle length   
854      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
855                                                         zwgt2    ! Secondary weights
856     
857      INTEGER ::  jic, jn, ji                      ! temporary integers
858      REAL(wp) :: za1, za2
859      !!----------------------------------------------------------------------
860
861      zwgt1(:) = 0._wp
862      zwgt2(:) = 0._wp
863
864      ! Set time index when averaged value is requested
865      IF (ll_fw) THEN
866         jic = nn_baro
867      ELSE
868         jic = 2 * nn_baro
869      ENDIF
870
871      ! Set primary weights:
872      IF (ll_av) THEN
873           ! Define simple boxcar window for primary weights
874           ! (width = nn_baro, centered around jic)     
875         SELECT CASE ( nn_bt_flt )
876              CASE( 0 )  ! No averaging
877                 zwgt1(jic) = 1._wp
878                 jpit = jic
879
880              CASE( 1 )  ! Boxcar, width = nn_baro
881                 DO jn = 1, 3*nn_baro
882                    za1 = ABS(float(jn-jic))/float(nn_baro) 
883                    IF (za1 < 0.5_wp) THEN
884                      zwgt1(jn) = 1._wp
885                      jpit = jn
886                    ENDIF
887                 ENDDO
888
889              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
890                 DO jn = 1, 3*nn_baro
891                    za1 = ABS(float(jn-jic))/float(nn_baro) 
892                    IF (za1 < 1._wp) THEN
893                      zwgt1(jn) = 1._wp
894                      jpit = jn
895                    ENDIF
896                 ENDDO
897              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
898         END SELECT
899
900      ELSE ! No time averaging
901         zwgt1(jic) = 1._wp
902         jpit = jic
903      ENDIF
904   
905      ! Set secondary weights
906      DO jn = 1, jpit
907        DO ji = jn, jpit
908             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
909        END DO
910      END DO
911
912      ! Normalize weigths:
913      za1 = 1._wp / SUM(zwgt1(1:jpit))
914      za2 = 1._wp / SUM(zwgt2(1:jpit))
915      DO jn = 1, jpit
916        zwgt1(jn) = zwgt1(jn) * za1
917        zwgt2(jn) = zwgt2(jn) * za2
918      END DO
919      !
920   END SUBROUTINE ts_wgt
921
922
923   SUBROUTINE ts_rst( kt, cdrw )
924      !!---------------------------------------------------------------------
925      !!                   ***  ROUTINE ts_rst  ***
926      !!
927      !! ** Purpose : Read or write time-splitting arrays in restart file
928      !!----------------------------------------------------------------------
929      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
930      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
931      !!----------------------------------------------------------------------
932      !
933      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
934         !                                   ! ---------------
935         IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN    !* Read the restart file
936            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
937            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
938            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
939            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
940            IF( .NOT.ln_bt_av ) THEN
941               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
942               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
943               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
944               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
945               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
946               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
947            ENDIF
948#if defined key_agrif
949            ! Read time integrated fluxes
950            IF ( .NOT.Agrif_Root() ) THEN
951               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
952               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
953            ENDIF
954#endif
955         ELSE                                   !* Start from rest
956            IF(lwp) WRITE(numout,*)
957            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
958            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
959            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
960            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
961#if defined key_agrif
962            IF ( .NOT.Agrif_Root() ) THEN
963               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
964            ENDIF
965#endif
966         ENDIF
967         !
968      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
969         !                                   ! -------------------
970         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
971         IF( lwxios ) CALL iom_swap(      cwxios_context          )
972         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
973         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
974         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
975         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
976         !
977         IF (.NOT.ln_bt_av) THEN
978            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
979            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
980            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
981            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
982            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
983            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
984         ENDIF
985#if defined key_agrif
986         ! Save time integrated fluxes
987         IF ( .NOT.Agrif_Root() ) THEN
988            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
989            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
990         ENDIF
991#endif
992         IF( lwxios ) CALL iom_swap(      cxios_context          )
993      ENDIF
994      !
995   END SUBROUTINE ts_rst
996
997
998   SUBROUTINE dyn_spg_ts_init
999      !!---------------------------------------------------------------------
1000      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1001      !!
1002      !! ** Purpose : Set time splitting options
1003      !!----------------------------------------------------------------------
1004      INTEGER  ::   ji ,jj              ! dummy loop indices
1005      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1006      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
1007      !!----------------------------------------------------------------------
1008      !
1009      ! Max courant number for ext. grav. waves
1010      !
1011      DO jj = 1, jpj
1012         DO ji =1, jpi
1013            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1014            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1015            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1016         END DO
1017      END DO
1018      !
1019      zcmax = MAXVAL( zcu(:,:) )
1020      CALL mpp_max( 'dynspg_ts', zcmax )
1021
1022      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1023      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1024     
1025      rdtbt = rdt / REAL( nn_baro , wp )
1026      zcmax = zcmax * rdtbt
1027      ! Print results
1028      IF(lwp) WRITE(numout,*)
1029      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1030      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1031      IF( ln_bt_auto ) THEN
1032         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
1033         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1034      ELSE
1035         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
1036      ENDIF
1037
1038      IF(ln_bt_av) THEN
1039         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
1040      ELSE
1041         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
1042      ENDIF
1043      !
1044      !
1045      IF(ln_bt_fw) THEN
1046         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1047      ELSE
1048         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1049      ENDIF
1050      !
1051#if defined key_agrif
1052      ! Restrict the use of Agrif to the forward case only
1053!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1054#endif
1055      !
1056      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1057      SELECT CASE ( nn_bt_flt )
1058         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1059         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1060         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1061         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1062      END SELECT
1063      !
1064      IF(lwp) WRITE(numout,*) ' '
1065      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1066      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1067      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1068      !
1069      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1070      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1071         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1072      ENDIF
1073      !
1074      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1075         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1076      ENDIF
1077      IF( zcmax>0.9_wp ) THEN
1078         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1079      ENDIF
1080      !
1081      !                             ! Allocate time-splitting arrays
1082      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1083      !
1084      !                             ! read restart when needed
1085      CALL ts_rst( nit000, 'READ' )
1086      !
1087      IF( lwxios ) THEN
1088! define variables in restart file when writing with XIOS
1089         CALL iom_set_rstw_var_active('ub2_b')
1090         CALL iom_set_rstw_var_active('vb2_b')
1091         CALL iom_set_rstw_var_active('un_bf')
1092         CALL iom_set_rstw_var_active('vn_bf')
1093         !
1094         IF (.NOT.ln_bt_av) THEN
1095            CALL iom_set_rstw_var_active('sshbb_e')
1096            CALL iom_set_rstw_var_active('ubb_e')
1097            CALL iom_set_rstw_var_active('vbb_e')
1098            CALL iom_set_rstw_var_active('sshb_e')
1099            CALL iom_set_rstw_var_active('ub_e')
1100            CALL iom_set_rstw_var_active('vb_e')
1101         ENDIF
1102#if defined key_agrif
1103         ! Save time integrated fluxes
1104         IF ( .NOT.Agrif_Root() ) THEN
1105            CALL iom_set_rstw_var_active('ub2_i_b')
1106            CALL iom_set_rstw_var_active('vb2_i_b')
1107         ENDIF
1108#endif
1109      ENDIF
1110      !
1111   END SUBROUTINE dyn_spg_ts_init
1112
1113   
1114   SUBROUTINE dyn_cor_2D_init( Kmm )
1115      !!---------------------------------------------------------------------
1116      !!                   ***  ROUTINE dyn_cor_2D_init  ***
1117      !!
1118      !! ** Purpose : Set time splitting options
1119      !! Set arrays to remove/compute coriolis trend.
1120      !! Do it once during initialization if volume is fixed, else at each long time step.
1121      !! Note that these arrays are also used during barotropic loop. These are however frozen
1122      !! although they should be updated in the variable volume case. Not a big approximation.
1123      !! To remove this approximation, copy lines below inside barotropic loop
1124      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1125      !!
1126      !! Compute zwz = f / ( height of the water colomn )
1127      !!----------------------------------------------------------------------
1128      INTEGER,  INTENT(in)         ::  Kmm  ! Time index
1129      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
1130      REAL(wp) ::   z1_ht
1131      REAL(wp), DIMENSION(jpi,jpj) :: zhf
1132      !!----------------------------------------------------------------------
1133      !
1134      SELECT CASE( nvor_scheme )
1135      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
1136         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
1137         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
1138            DO jj = 1, jpjm1
1139               DO ji = 1, jpim1
1140                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
1141                       &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp 
1142                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1143               END DO
1144            END DO
1145         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
1146            DO jj = 1, jpjm1
1147               DO ji = 1, jpim1
1148                  zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
1149                       &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
1150                       &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
1151                       &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
1152                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1153               END DO
1154            END DO
1155         END SELECT
1156         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1157         !
1158         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1159         DO jj = 2, jpj
1160            DO ji = 2, jpi
1161               ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1162               ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1163               ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1164               ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1165            END DO
1166         END DO
1167         !
1168      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
1169         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
1170         DO jj = 2, jpj
1171            DO ji = 2, jpi
1172               z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1173               ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
1174               ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
1175               ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1176               ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
1177            END DO
1178         END DO
1179         !
1180      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
1181         !
1182         zwz(:,:) = 0._wp
1183         zhf(:,:) = 0._wp
1184         
1185         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1186!!gm    A priori a better value should be something like :
1187!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1188!!gm                     divided by the sum of the corresponding mask
1189!!gm
1190!!           
1191         IF( .NOT.ln_sco ) THEN
1192 
1193   !!gm  agree the JC comment  : this should be done in a much clear way
1194 
1195   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1196   !     Set it to zero for the time being
1197   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
1198   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
1199   !              ENDIF
1200   !              zhf(:,:) = gdepw_0(:,:,jk+1)
1201            !
1202         ELSE
1203            !
1204            !zhf(:,:) = hbatf(:,:)
1205            DO jj = 1, jpjm1
1206               DO ji = 1, jpim1
1207                  zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
1208                       &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
1209                       &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
1210                       &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
1211               END DO
1212            END DO
1213         ENDIF
1214         !
1215         DO jj = 1, jpjm1
1216            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1217         END DO
1218         !
1219         DO jk = 1, jpkm1
1220            DO jj = 1, jpjm1
1221               zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1222            END DO
1223         END DO
1224         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1225         ! JC: TBC. hf should be greater than 0
1226         DO jj = 1, jpj
1227            DO ji = 1, jpi
1228               IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
1229            END DO
1230         END DO
1231         zwz(:,:) = ff_f(:,:) * zwz(:,:)
1232      END SELECT
1233     
1234   END SUBROUTINE dyn_cor_2d_init
1235
1236
1237
1238   SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
1239      !!---------------------------------------------------------------------
1240      !!                   ***  ROUTINE dyn_cor_2d  ***
1241      !!
1242      !! ** Purpose : Compute u and v coriolis trends
1243      !!----------------------------------------------------------------------
1244      INTEGER  ::   ji ,jj                             ! dummy loop indices
1245      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
1246      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phu, phv, punb, pvnb, zhU, zhV
1247      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
1248      !!----------------------------------------------------------------------
1249      SELECT CASE( nvor_scheme )
1250      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
1251         DO jj = 2, jpjm1
1252            DO ji = 2, jpim1
1253               z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1254               z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1255               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
1256                  &               * (  e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
1257                  &                  + e1e2t(ji  ,jj)*ht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
1258                  !
1259               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1260                  &               * (  e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
1261                  &                  + e1e2t(ji,jj  )*ht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
1262            END DO 
1263         END DO 
1264         !         
1265      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
1266         DO jj = 2, jpjm1
1267            DO ji = fs_2, fs_jpim1   ! vector opt.
1268               zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1269               zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1270               zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1271               zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1272               ! energy conserving formulation for planetary vorticity term
1273               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1274               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1275            END DO
1276         END DO
1277         !
1278      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
1279         DO jj = 2, jpjm1
1280            DO ji = fs_2, fs_jpim1   ! vector opt.
1281               zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
1282                 &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1283               zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
1284                 &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1285               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1286               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1287            END DO
1288         END DO
1289         !
1290      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
1291         DO jj = 2, jpjm1
1292            DO ji = fs_2, fs_jpim1   ! vector opt.
1293               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
1294                &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
1295                &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
1296                &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1297               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1298                &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
1299                &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
1300                &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
1301            END DO
1302         END DO
1303         !
1304      END SELECT
1305      !
1306   END SUBROUTINE dyn_cor_2D
1307
1308
1309   SUBROUTINE wad_tmsk( pssh, ptmsk )
1310      !!----------------------------------------------------------------------
1311      !!                  ***  ROUTINE wad_lmt  ***
1312      !!                   
1313      !! ** Purpose :   set wetting & drying mask at tracer points
1314      !!              for the current barotropic sub-step
1315      !!
1316      !! ** Method  :   ???
1317      !!
1318      !! ** Action  :  ptmsk : wetting & drying t-mask
1319      !!----------------------------------------------------------------------
1320      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
1321      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
1322      !
1323      INTEGER  ::   ji, jj   ! dummy loop indices
1324      !!----------------------------------------------------------------------
1325      !
1326      IF( ln_wd_dl_rmp ) THEN     
1327         DO jj = 1, jpj
1328            DO ji = 1, jpi                   
1329               IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
1330                  !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
1331                  ptmsk(ji,jj) = 1._wp
1332               ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
1333                  ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
1334               ELSE
1335                  ptmsk(ji,jj) = 0._wp
1336               ENDIF
1337            END DO
1338         END DO
1339      ELSE 
1340         DO jj = 1, jpj
1341            DO ji = 1, jpi                             
1342               IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
1343               ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
1344               ENDIF
1345            END DO
1346         END DO
1347      ENDIF
1348      !
1349   END SUBROUTINE wad_tmsk
1350
1351
1352   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1353      !!----------------------------------------------------------------------
1354      !!                  ***  ROUTINE wad_lmt  ***
1355      !!                   
1356      !! ** Purpose :   set wetting & drying mask at tracer points
1357      !!              for the current barotropic sub-step
1358      !!
1359      !! ** Method  :   ???
1360      !!
1361      !! ** Action  :  ptmsk : wetting & drying t-mask
1362      !!----------------------------------------------------------------------
1363      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
1364      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
1365      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
1366      !
1367      INTEGER  ::   ji, jj   ! dummy loop indices
1368      !!----------------------------------------------------------------------
1369      !
1370      DO jj = 1, jpj
1371         DO ji = 1, jpim1   ! not jpi-column
1372            IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
1373            ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj) 
1374            ENDIF
1375            phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1376            pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1377         END DO
1378      END DO
1379      !
1380      DO jj = 1, jpjm1   ! not jpj-row
1381         DO ji = 1, jpi
1382            IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
1383            ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1) 
1384            ENDIF
1385            phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
1386            pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1387         END DO
1388      END DO
1389      !
1390   END SUBROUTINE wad_Umsk
1391
1392
1393   SUBROUTINE wad_spg( pshn, zcpx, zcpy )
1394      !!---------------------------------------------------------------------
1395      !!                   ***  ROUTINE  wad_sp  ***
1396      !!
1397      !! ** Purpose :
1398      !!----------------------------------------------------------------------
1399      INTEGER  ::   ji ,jj               ! dummy loop indices
1400      LOGICAL  ::   ll_tmp1, ll_tmp2
1401      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
1402      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1403      !!----------------------------------------------------------------------
1404      DO jj = 2, jpjm1
1405         DO ji = 2, jpim1 
1406            ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
1407                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
1408                 &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji+1,jj) + ht_0(ji+1,jj) )  &
1409                 &                                                         > rn_wdmin1 + rn_wdmin2
1410            ll_tmp2 = ( ABS( pshn(ji+1,jj)            -  pshn(ji  ,jj))  > 1.E-12 ).AND.( &
1411                 &      MAX(   pshn(ji,jj)              ,  pshn(ji+1,jj) ) >                &
1412                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1413            IF(ll_tmp1) THEN
1414               zcpx(ji,jj) = 1.0_wp
1415            ELSEIF(ll_tmp2) THEN
1416               ! no worries about  pshn(ji+1,jj) -  pshn(ji  ,jj) = 0, it won't happen ! here
1417               zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1418                    &           / (pshn(ji+1,jj) - pshn(ji  ,jj)) )
1419               zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1420            ELSE
1421               zcpx(ji,jj) = 0._wp
1422            ENDIF
1423            !
1424            ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji,jj+1) ) >                &
1425                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
1426                 &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji,jj+1) + ht_0(ji,jj+1) )  &
1427                 &                                                       > rn_wdmin1 + rn_wdmin2
1428            ll_tmp2 = ( ABS( pshn(ji,jj)              -  pshn(ji,jj+1))  > 1.E-12 ).AND.( &
1429                 &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
1430                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1431           
1432            IF(ll_tmp1) THEN
1433               zcpy(ji,jj) = 1.0_wp
1434            ELSE IF(ll_tmp2) THEN
1435               ! no worries about  pshn(ji,jj+1) -  pshn(ji,jj  ) = 0, it won't happen ! here
1436               zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1437                    &           / (pshn(ji,jj+1) - pshn(ji,jj  )) )
1438               zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
1439            ELSE
1440               zcpy(ji,jj) = 0._wp
1441            ENDIF
1442         END DO
1443      END DO
1444           
1445   END SUBROUTINE wad_spg
1446     
1447
1448
1449   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
1450      !!----------------------------------------------------------------------
1451      !!                  ***  ROUTINE dyn_drg_init  ***
1452      !!                   
1453      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
1454      !!              the baroclinic part of the barotropic RHS
1455      !!              - compute the barotropic drag coefficients
1456      !!
1457      !! ** Method  :   computation done over the INNER domain only
1458      !!----------------------------------------------------------------------
1459      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
1460      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
1461      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
1462      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
1463      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients
1464      !
1465      INTEGER  ::   ji, jj   ! dummy loop indices
1466      INTEGER  ::   ikbu, ikbv, iktu, iktv
1467      REAL(wp) ::   zztmp
1468      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
1469      !!----------------------------------------------------------------------
1470      !
1471      !                    !==  Set the barotropic drag coef.  ==!
1472      !
1473      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
1474         
1475         DO jj = 2, jpjm1
1476            DO ji = 2, jpim1     ! INNER domain
1477               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1478               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1479            END DO
1480         END DO
1481      ELSE                          ! bottom friction only
1482         DO jj = 2, jpjm1
1483            DO ji = 2, jpim1  ! INNER domain
1484               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1485               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1486            END DO
1487         END DO
1488      ENDIF
1489      !
1490      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
1491      !
1492      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
1493         
1494         DO jj = 2, jpjm1
1495            DO ji = 2, jpim1  ! INNER domain
1496               ikbu = mbku(ji,jj)       
1497               ikbv = mbkv(ji,jj)   
1498               zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1499               zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1500            END DO
1501         END DO
1502      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
1503         
1504         DO jj = 2, jpjm1
1505            DO ji = 2, jpim1   ! INNER domain
1506               ikbu = mbku(ji,jj)       
1507               ikbv = mbkv(ji,jj)   
1508               zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1509               zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1510            END DO
1511         END DO
1512      ENDIF
1513      !
1514      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
1515         zztmp = -1._wp / rdtbt
1516         DO jj = 2, jpjm1
1517            DO ji = 2, jpim1    ! INNER domain
1518               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
1519                    &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
1520               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
1521                    &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
1522            END DO
1523         END DO
1524      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1525         
1526         DO jj = 2, jpjm1
1527            DO ji = 2, jpim1    ! INNER domain
1528               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)
1529               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)
1530            END DO
1531         END DO
1532      END IF
1533      !
1534      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
1535      !
1536      IF( ln_isfcav ) THEN
1537         !
1538         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
1539           
1540            DO jj = 2, jpjm1
1541               DO ji = 2, jpim1   ! INNER domain
1542                  iktu = miku(ji,jj)
1543                  iktv = mikv(ji,jj)
1544                  zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1545                  zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1546               END DO
1547            END DO
1548         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
1549           
1550            DO jj = 2, jpjm1
1551               DO ji = 2, jpim1      ! INNER domain
1552                  iktu = miku(ji,jj)
1553                  iktv = mikv(ji,jj)
1554                  zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1555                  zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1556               END DO
1557            END DO
1558         ENDIF
1559         !
1560         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1561         
1562         DO jj = 2, jpjm1
1563            DO ji = 2, jpim1    ! INNER domain
1564               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)
1565               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)
1566            END DO
1567         END DO
1568         !
1569      ENDIF
1570      !
1571   END SUBROUTINE dyn_drg_init
1572
1573   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
1574      &                      za0, za1, za2, za3 )   ! ==> out
1575      !!----------------------------------------------------------------------
1576      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
1577      LOGICAL ,INTENT(in   ) ::   ll_init              !
1578      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient
1579      !
1580      REAL(wp) ::   zepsilon, zgamma                   !   -      -
1581      !!----------------------------------------------------------------------
1582      !                             ! set Half-step back interpolation coefficient
1583      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
1584         za0 = 1._wp                       
1585         za1 = 0._wp                           
1586         za2 = 0._wp
1587         za3 = 0._wp
1588      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1589         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
1590         za1 =-0.1666666666666_wp                 ! za1 = gam
1591         za2 = 0.0833333333333_wp                 ! za2 = eps
1592         za3 = 0._wp             
1593      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1594         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion 
1595            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
1596            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
1597            za2 = 0.088_wp                        ! za2 = gam
1598            za3 = 0.013_wp                        ! za3 = eps
1599         ELSE                                 ! no time diffusion
1600            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1601            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1602            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1603            za1 = 1._wp - za0 - zgamma - zepsilon
1604            za2 = zgamma
1605            za3 = zepsilon
1606         ENDIF
1607      ENDIF
1608   END SUBROUTINE ts_bck_interp
1609
1610
1611   !!======================================================================
1612END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.