source: NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DYN/dynspg_ts.F90 @ 12117

Last change on this file since 12117 was 12117, checked in by smueller, 9 months ago

Merging of the developments in /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides@12096 with respect to /NEMO/trunk@12072 into /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis (tickets #2175 and #2194)

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