source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DYN/dynspg_ts.F90 @ 12147

Last change on this file since 12147 was 12147, checked in by gsamson, 10 months ago

update dev_ASINTER-01-05_merged to trunk@r12072 + bugfix for ICE_AGRIF test case

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