source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 12616

Last change on this file since 12616 was 12616, checked in by techene, 11 months ago

all: add e3 substitute (sometimes it requires to add ze3t/u/v/w) and limit precompiled files lines to about 130 character, OCE/ASM/asminc.F90, OCE/DOM/domzgr_substitute.h90, OCE/ISF/isfcpl.F90, OCE/SBC/sbcice_cice, OCE/CRS/crsini.F90 : add key_LF

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