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

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

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90 @ 13220

Last change on this file since 13220 was 13220, checked in by orioltp, 4 years ago

dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: updating from trunk r13218

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