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/MOI/NEMO_4.03_IODRAG/src/OCE/DYN – NEMO

source: NEMO/branches/MOI/NEMO_4.03_IODRAG/src/OCE/DYN/dynspg_ts.F90 @ 12953

Last change on this file since 12953 was 12953, checked in by jchanut, 4 years ago

Add ice-ocean implicit drag as namelist option ln_drgice_imp, #2468

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