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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

  • Property svn:keywords set to Id
File size: 63.6 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_ts  ***
4   !! Ocean dynamics:  surface pressure gradient trend, split-explicit scheme
5   !!======================================================================
6   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
7   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
8   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
9   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
10   !!              -   ! 2008-01  (R. Benshila)  change averaging method
11   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
12   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
13   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
14   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
15   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
16   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
17   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
18   !!---------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
22   !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
23   !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
24   !!   ts_rst         : read/write time-splitting fields in restart file
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce         ! surface boundary condition: ocean
29   USE zdf_oce         ! Bottom friction coefts
30   USE sbcisf          ! ice shelf variable (fwfisf)
31   USE sbcapr          ! surface boundary condition: atmospheric pressure
32   USE dynadv    , ONLY: ln_dynadv_vec
33   USE phycst          ! physical constants
34   USE dynvor          ! vorticity term
35   USE wet_dry         ! wetting/drying flux limter
36   USE bdy_oce         ! open boundary
37   USE bdytides        ! open boundary condition data
38   USE bdydyn2d        ! open boundary conditions on barotropic variables
39   USE sbctide         ! tides
40   USE updtide         ! tide potential
41   USE sbcwave         ! surface wave
42   !
43   USE in_out_manager  ! I/O manager
44   USE lib_mpp         ! distributed memory computing library
45   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
46   USE prtctl          ! Print control
47   USE iom             ! IOM library
48   USE restart         ! only for lrst_oce
49   USE wrk_nemo        ! Memory Allocation
50   USE timing          ! Timing   
51   USE diatmb          ! Top,middle,bottom output
52#if defined key_agrif
53   USE agrif_opa_interp ! agrif
54#endif
55#if defined key_asminc   
56   USE asminc          ! Assimilation increment
57#endif
58
59
60   IMPLICIT NONE
61   PRIVATE
62
63   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90
64   PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
65   PUBLIC dyn_spg_ts_init   !    "      "     "    "
66   PUBLIC ts_rst            !    "      "     "    "
67
68   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
69   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
70
71   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   wgtbtp1, wgtbtp2   !: 1st & 2nd weights used in time filtering of barotropic fields
72
73   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff_f/h at F points
74   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   !: triad of coriolis parameter
75   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   !: (only used with een vorticity scheme)
76
77   !! Time filtered arrays at baroclinic time step:
78   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv     !: Advection vel. at "now" barocl. step
79
80   !! * Substitutions
81#  include "vectopt_loop_substitute.h90"
82   !!----------------------------------------------------------------------
83   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
84   !! $Id$
85   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
86   !!----------------------------------------------------------------------
87CONTAINS
88
89   INTEGER FUNCTION dyn_spg_ts_alloc()
90      !!----------------------------------------------------------------------
91      !!                  ***  routine dyn_spg_ts_alloc  ***
92      !!----------------------------------------------------------------------
93      INTEGER :: ierr(3)
94      !!----------------------------------------------------------------------
95      ierr(:) = 0
96      !
97      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
98      !
99      IF( ln_dynvor_een )   ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
100         &                            ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
101         !
102      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
103      !
104      dyn_spg_ts_alloc = MAXVAL( ierr(:) )
105      !
106      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
107      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
108      !
109   END FUNCTION dyn_spg_ts_alloc
110
111
112   SUBROUTINE dyn_spg_ts( kt )
113      !!----------------------------------------------------------------------
114      !!
115      !! ** Purpose : - Compute the now trend due to the explicit time stepping
116      !!              of the quasi-linear barotropic system, and add it to the
117      !!              general momentum trend.
118      !!
119      !! ** Method  : - split-explicit schem (time splitting) :
120      !!      Barotropic variables are advanced from internal time steps
121      !!      "n"   to "n+1" if ln_bt_fw=T
122      !!      or from
123      !!      "n-1" to "n+1" if ln_bt_fw=F
124      !!      thanks to a generalized forward-backward time stepping (see ref. below).
125      !!
126      !! ** Action :
127      !!      -Update the filtered free surface at step "n+1"      : ssha
128      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
129      !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
130      !!      These are used to advect tracers and are compliant with discrete
131      !!      continuity equation taken at the baroclinic time steps. This
132      !!      ensures tracers conservation.
133      !!      - (ua, va) momentum trend updated with barotropic component.
134      !!
135      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
136      !!---------------------------------------------------------------------
137      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
138      !
139      LOGICAL  ::   ll_fw_start        ! if true, forward integration
140      LOGICAL  ::   ll_init             ! if true, special startup of 2d equations
141      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D
142      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
143      INTEGER  ::   ikbu, ikbv, noffset      ! local integers
144      INTEGER  ::   iktu, iktv               ! local integers
145      REAL(wp) ::   zmdi
146      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
147      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      -
148      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      -
149      REAL(wp) ::   zu_spg, zv_spg              !   -      -
150      REAL(wp) ::   zhura, zhvra          !   -      -
151      REAL(wp) ::   za0, za1, za2, za3    !   -      -
152      !
153      REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e
154      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
155      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv
156      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
157      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
158      REAL(wp), POINTER, DIMENSION(:,:) :: zhf
159      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef.
160      !!----------------------------------------------------------------------
161      !
162      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
163      !
164      !                                         !* Allocate temporary arrays
165      CALL wrk_alloc( jpi,jpj,   zsshp2_e, zhdiv )
166      CALL wrk_alloc( jpi,jpj,   zu_trd, zv_trd)
167      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc)
168      CALL wrk_alloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e)
169      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  )
170      CALL wrk_alloc( jpi,jpj,   zhf )
171      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
172      !
173      zmdi=1.e+20                               !  missing data indicator for masking
174      !                                         !* Local constant initialization
175      z1_12 = 1._wp / 12._wp 
176      z1_8  = 0.125_wp                                   
177      z1_4  = 0.25_wp
178      z1_2  = 0.5_wp     
179      zraur = 1._wp / rau0
180      !                                            ! reciprocal of baroclinic time step
181      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
182      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt
183      ENDIF
184      z1_2dt_b = 1.0_wp / z2dt_bf 
185      !
186      ll_init     = ln_bt_av                       ! if no time averaging, then no specific restart
187      ll_fw_start = .FALSE.
188      !                                            ! time offset in steps for bdy data update
189      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro
190      ELSE                       ;   noffset =   0 
191      ENDIF
192      !
193      IF( kt == nit000 ) THEN                !* initialisation
194         !
195         IF(lwp) WRITE(numout,*)
196         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
197         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
198         IF(lwp) WRITE(numout,*)
199         !
200         IF( neuler == 0 )   ll_init=.TRUE.
201         !
202         IF( ln_bt_fw .OR. neuler == 0 ) THEN
203            ll_fw_start =.TRUE.
204            noffset     = 0
205         ELSE
206            ll_fw_start =.FALSE.
207         ENDIF
208         !
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      ! Set arrays to remove/compute coriolis trend.
215      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
216      ! Note that these arrays are also used during barotropic loop. These are however frozen
217      ! although they should be updated in the variable volume case. Not a big approximation.
218      ! To remove this approximation, copy lines below inside barotropic loop
219      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
220      !
221      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
222         IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==!
223            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
224            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
225               DO jj = 1, jpjm1
226                  DO ji = 1, jpim1
227                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
228                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp 
229                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
230                  END DO
231               END DO
232            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
233               DO jj = 1, jpjm1
234                  DO ji = 1, jpim1
235                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     &
236                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   &
237                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
238                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
239                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
240                  END DO
241               END DO
242            END SELECT
243            CALL lbc_lnk( zwz, 'F', 1._wp )
244            !
245            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
246            DO jj = 2, jpj
247               DO ji = 2, jpi
248                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
249                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
250                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
251                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
252               END DO
253            END DO
254            !
255         ELSE                                !== all other schemes (ENE, ENS, MIX)
256            zwz(:,:) = 0._wp
257            zhf(:,:) = 0._wp
258           
259!!gm  assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed
260!!gm    A priori a better value should be something like :
261!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
262!!gm                     divided by the sum of the corresponding mask
263!!gm
264!!           
265              IF ( .not. ln_sco ) THEN
266 
267   !!gm  agree the JC comment  : this should be done in a much clear way
268 
269   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
270   !     Set it to zero for the time being
271   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
272   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
273   !              ENDIF
274   !              zhf(:,:) = gdepw_0(:,:,jk+1)
275               ELSE
276                 !zhf(:,:) = hbatf(:,:)
277                 DO jj = 1, jpjm1
278                   DO ji = 1, jpim1
279                     zhf(ji,jj) = MAX( 0._wp,                                &
280                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     &
281                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     &
282                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     &
283                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   &
284                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +&
285                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +&
286                                &   rsmall  ) )
287                   END DO
288                 END DO
289              END IF
290 
291              DO jj = 1, jpjm1
292                 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
293              END DO
294!!gm end
295
296            DO jk = 1, jpkm1
297               DO jj = 1, jpjm1
298                  zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
299               END DO
300            END DO
301            CALL lbc_lnk( zhf, 'F', 1._wp )
302            ! JC: TBC. hf should be greater than 0
303            DO jj = 1, jpj
304               DO ji = 1, jpi
305                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
306               END DO
307            END DO
308            zwz(:,:) = ff_f(:,:) * zwz(:,:)
309         ENDIF
310      ENDIF
311      !
312      ! If forward start at previous time step, and centered integration,
313      ! then update averaging weights:
314      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
315         ll_fw_start=.FALSE.
316         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
317      ENDIF
318                         
319      ! -----------------------------------------------------------------------------
320      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
321      ! -----------------------------------------------------------------------------
322      !     
323      !
324      !                                   !* e3*d/dt(Ua) (Vertically integrated)
325      !                                   ! --------------------------------------------------
326      zu_frc(:,:) = 0._wp
327      zv_frc(:,:) = 0._wp
328      !
329      DO jk = 1, jpkm1
330         zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
331         zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
332      END DO
333      !
334      zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
335      zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
336      !
337      !
338      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
339      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
340         DO jj = 2, jpjm1
341            DO ji = fs_2, fs_jpim1   ! vector opt.
342               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
343               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
344            END DO
345         END DO
346      END DO
347     
348!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
349!!gm  Is it correct to do so ?   I think so...
350     
351     
352      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
353      !                                   ! --------------------------------------------------------
354      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes
355      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
356      !
357      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
358         DO jj = 2, jpjm1
359            DO ji = fs_2, fs_jpim1   ! vector opt.
360               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
361               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
362               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
363               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
364               ! energy conserving formulation for planetary vorticity term
365               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
366               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
367            END DO
368         END DO
369         !
370      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
371         DO jj = 2, jpjm1
372            DO ji = fs_2, fs_jpim1   ! vector opt.
373               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
374                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
375               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
376                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
377               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
378               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
379            END DO
380         END DO
381         !
382      ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme
383         DO jj = 2, jpjm1
384            DO ji = fs_2, fs_jpim1   ! vector opt.
385               zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
386                &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
387                &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
388                &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
389               zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
390                &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
391                &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
392                &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
393            END DO
394         END DO
395         !
396      ENDIF 
397      !
398      !                                   !* Right-Hand-Side of the barotropic momentum equation
399      !                                   ! ----------------------------------------------------
400      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
401        IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters
402           DO jj = 2, jpjm1
403              DO ji = 2, jpim1 
404                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
405                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
406                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
407                     &                                                         > rn_wdmin1 + rn_wdmin2
408                ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( &
409                     &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
410                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
411   
412                IF(ll_tmp1) THEN
413                  zcpx(ji,jj) = 1.0_wp
414                ELSE IF(ll_tmp2) THEN
415                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
416                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
417                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
418                ELSE
419                  zcpx(ji,jj) = 0._wp
420                END IF
421         
422                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
423                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
424                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
425                     &                                                         > rn_wdmin1 + rn_wdmin2
426                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( &
427                     &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
428                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
429   
430                IF(ll_tmp1) THEN
431                  zcpy(ji,jj) = 1.0_wp
432                ELSE IF(ll_tmp2) THEN
433                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
434                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
435                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
436                ELSE
437                  zcpy(ji,jj) = 0._wp
438                END IF
439              END DO
440           END DO
441 
442           DO jj = 2, jpjm1
443              DO ji = 2, jpim1
444                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
445                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj)
446                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
447                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj)
448              END DO
449           END DO
450
451         ELSE
452
453           DO jj = 2, jpjm1
454              DO ji = fs_2, fs_jpim1   ! vector opt.
455                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
456                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
457              END DO
458           END DO
459        ENDIF
460
461      ENDIF
462
463      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
464         DO ji = fs_2, fs_jpim1
465             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
466             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
467          END DO
468      END DO 
469      !
470      !                 ! Add bottom stress contribution from baroclinic velocities:     
471      IF (ln_bt_fw) THEN
472         DO jj = 2, jpjm1                         
473            DO ji = fs_2, fs_jpim1   ! vector opt.
474               ikbu = mbku(ji,jj)       
475               ikbv = mbkv(ji,jj)   
476               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
477               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
478            END DO
479         END DO
480      ELSE
481         DO jj = 2, jpjm1
482            DO ji = fs_2, fs_jpim1   ! vector opt.
483               ikbu = mbku(ji,jj)       
484               ikbv = mbkv(ji,jj)   
485               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
486               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
487            END DO
488         END DO
489      ENDIF
490      !
491      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
492      IF( ln_wd ) THEN
493        zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:)
494        zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:)
495      ELSE
496        zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:)
497        zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:)
498      END IF
499      !
500      !                                         ! Add top stress contribution from baroclinic velocities:     
501      IF( ln_bt_fw ) THEN
502         DO jj = 2, jpjm1
503            DO ji = fs_2, fs_jpim1   ! vector opt.
504               iktu = miku(ji,jj)
505               iktv = mikv(ji,jj)
506               zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
507               zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
508            END DO
509         END DO
510      ELSE
511         DO jj = 2, jpjm1
512            DO ji = fs_2, fs_jpim1   ! vector opt.
513               iktu = miku(ji,jj)
514               iktv = mikv(ji,jj)
515               zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
516               zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
517            END DO
518         END DO
519      ENDIF
520      !
521      ! Note that the "unclipped" top friction parameter is used even with explicit drag
522      zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:)
523      zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:)
524      !       
525      IF (ln_bt_fw) THEN                        ! Add wind forcing
526         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:)
527         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:)
528      ELSE
529         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:)
530         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:)
531      ENDIF 
532      !
533      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
534         IF (ln_bt_fw) THEN
535            DO jj = 2, jpjm1             
536               DO ji = fs_2, fs_jpim1   ! vector opt.
537                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
538                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
539                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
540                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
541               END DO
542            END DO
543         ELSE
544            DO jj = 2, jpjm1             
545               DO ji = fs_2, fs_jpim1   ! vector opt.
546                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
547                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
548                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
549                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
550                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
551                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
552               END DO
553            END DO
554         ENDIF
555      ENDIF
556      !                                   !* Right-Hand-Side of the barotropic ssh equation
557      !                                   ! -----------------------------------------------
558      !                                         ! Surface net water flux and rivers
559      IF (ln_bt_fw) THEN
560         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
561      ELSE
562         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
563                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
564      ENDIF
565      !
566      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary
567         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
568      ENDIF
569      !
570#if defined key_asminc
571      !                                         ! Include the IAU weighted SSH increment
572      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
573         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
574      ENDIF
575#endif
576      !                                   !* Fill boundary data arrays for AGRIF
577      !                                   ! ------------------------------------
578#if defined key_agrif
579         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
580#endif
581      !
582      ! -----------------------------------------------------------------------
583      !  Phase 2 : Integration of the barotropic equations
584      ! -----------------------------------------------------------------------
585      !
586      !                                             ! ==================== !
587      !                                             !    Initialisations   !
588      !                                             ! ==================== ! 
589      ! Initialize barotropic variables:     
590      IF( ll_init )THEN
591         sshbb_e(:,:) = 0._wp
592         ubb_e  (:,:) = 0._wp
593         vbb_e  (:,:) = 0._wp
594         sshb_e (:,:) = 0._wp
595         ub_e   (:,:) = 0._wp
596         vb_e   (:,:) = 0._wp
597      ENDIF
598
599      !
600      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
601         sshn_e(:,:) =    sshn(:,:)           
602         un_e  (:,:) =    un_b(:,:)           
603         vn_e  (:,:) =    vn_b(:,:)
604         !
605         hu_e  (:,:) =    hu_n(:,:)       
606         hv_e  (:,:) =    hv_n(:,:) 
607         hur_e (:,:) = r1_hu_n(:,:)   
608         hvr_e (:,:) = r1_hv_n(:,:)
609      ELSE                                ! CENTRED integration: start from BEFORE fields
610         sshn_e(:,:) =    sshb(:,:)
611         un_e  (:,:) =    ub_b(:,:)         
612         vn_e  (:,:) =    vb_b(:,:)
613         !
614         hu_e  (:,:) =    hu_b(:,:)       
615         hv_e  (:,:) =    hv_b(:,:) 
616         hur_e (:,:) = r1_hu_b(:,:)   
617         hvr_e (:,:) = r1_hv_b(:,:)
618      ENDIF
619      !
620      !
621      !
622      ! Initialize sums:
623      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
624      va_b  (:,:) = 0._wp
625      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
626      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
627      vn_adv(:,:) = 0._wp
628      !                                             ! ==================== !
629      DO jn = 1, icycle                             !  sub-time-step loop  !
630         !                                          ! ==================== !
631         !                                                !* Update the forcing (BDY and tides)
632         !                                                !  ------------------
633         ! Update only tidal forcing at open boundaries
634#if defined key_tide
635         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
636         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   )
637#endif
638         !
639         ! Set extrapolation coefficients for predictor step:
640         IF ((jn<3).AND.ll_init) THEN      ! Forward           
641           za1 = 1._wp                                         
642           za2 = 0._wp                       
643           za3 = 0._wp                       
644         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
645           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
646           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
647           za3 =  0.281105_wp              ! za3 = bet
648         ENDIF
649
650         ! Extrapolate barotropic velocities at step jit+0.5:
651         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
652         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
653
654         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
655            !                                             !  ------------------
656            ! Extrapolate Sea Level at step jit+0.5:
657            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
658            !
659            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
660               DO ji = 2, fs_jpim1   ! Vector opt.
661                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
662                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
663                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
664                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
665                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
666                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
667               END DO
668            END DO
669            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
670            !
671            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
672            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
673         ELSE
674            zhup2_e (:,:) = hu_n(:,:)
675            zhvp2_e (:,:) = hv_n(:,:)
676         ENDIF
677         !                                                !* after ssh
678         !                                                !  -----------
679         ! One should enforce volume conservation at open boundaries here
680         ! considering fluxes below:
681         !
682         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
683         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
684         !
685#if defined key_agrif
686         ! Set fluxes during predictor step to ensure volume conservation
687         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
688            IF((nbondi == -1).OR.(nbondi == 2)) THEN
689               DO jj=1,jpj
690                  zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
691               END DO
692            ENDIF
693            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
694               DO jj=1,jpj
695                  zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
696               END DO
697            ENDIF
698            IF((nbondj == -1).OR.(nbondj == 2)) THEN
699               DO ji=1,jpi
700                  zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
701               END DO
702            ENDIF
703            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
704               DO ji=1,jpi
705                  zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
706               END DO
707            ENDIF
708         ENDIF
709#endif
710         IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
711         !
712         ! Sum over sub-time-steps to compute advective velocities
713         za2 = wgtbtp2(jn)
714         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
715         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
716         !
717         ! Set next sea level:
718         DO jj = 2, jpjm1                                 
719            DO ji = fs_2, fs_jpim1   ! vector opt.
720               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
721                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
722            END DO
723         END DO
724         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
725         
726         CALL lbc_lnk( ssha_e, 'T',  1._wp )
727
728         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
729         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
730#if defined key_agrif
731         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
732#endif
733         
734         ! Sea Surface Height at u-,v-points (vvl case only)
735         IF( .NOT.ln_linssh ) THEN                               
736            DO jj = 2, jpjm1
737               DO ji = 2, jpim1      ! NO Vector Opt.
738                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
739                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
740                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
741                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
742                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
743                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
744               END DO
745            END DO
746            CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
747         ENDIF   
748         !                                 
749         ! Half-step back interpolation of SSH for surface pressure computation:
750         !----------------------------------------------------------------------
751         IF ((jn==1).AND.ll_init) THEN
752           za0=1._wp                        ! Forward-backward
753           za1=0._wp                           
754           za2=0._wp
755           za3=0._wp
756         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
757           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
758           za1=-0.1666666666666_wp          ! za1 = gam
759           za2= 0.0833333333333_wp          ! za2 = eps
760           za3= 0._wp             
761         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
762           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
763           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
764           za2=0.088_wp                     ! za2 = gam
765           za3=0.013_wp                     ! za3 = eps
766         ENDIF
767         !
768         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
769          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
770         IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters
771           DO jj = 2, jpjm1
772              DO ji = 2, jpim1 
773                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
774                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) .AND.            &
775                     &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &
776                     &                                                             > rn_wdmin1 + rn_wdmin2
777                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
778                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
779                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
780   
781                IF(ll_tmp1) THEN
782                  zcpx(ji,jj) = 1.0_wp
783                ELSE IF(ll_tmp2) THEN
784                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here
785                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +    ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &
786                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) )
787                ELSE
788                  zcpx(ji,jj) = 0._wp
789                END IF
790         
791                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
792                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) .AND.            &
793                     &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &
794                     &                                                             > rn_wdmin1 + rn_wdmin2
795                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
796                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
797                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
798   
799                IF(ll_tmp1) THEN
800                  zcpy(ji,jj) = 1.0_wp
801                ELSE IF(ll_tmp2) THEN
802                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
803                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +    ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &
804                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) )
805                ELSE
806                  zcpy(ji,jj) = 0._wp
807                END IF
808              END DO
809           END DO
810         END IF
811         !
812         ! Compute associated depths at U and V points:
813         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN   !* Vector form
814            !                                       
815            DO jj = 2, jpjm1                           
816               DO ji = 2, jpim1
817                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
818                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
819                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
820                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
821                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
822                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
823                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
824                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
825               END DO
826            END DO
827
828         ENDIF
829         !
830         ! Add Coriolis trend:
831         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
832         ! at each time step. We however keep them constant here for optimization.
833         ! Recall that zwx and zwy arrays hold fluxes at this stage:
834         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
835         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
836         !
837         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN     !==  energy conserving or mixed scheme  ==!
838            DO jj = 2, jpjm1
839               DO ji = fs_2, fs_jpim1   ! vector opt.
840                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
841                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
842                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
843                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
844                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
845                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
846               END DO
847            END DO
848            !
849         ELSEIF ( ln_dynvor_ens ) THEN                   !==  enstrophy conserving scheme  ==!
850            DO jj = 2, jpjm1
851               DO ji = fs_2, fs_jpim1   ! vector opt.
852                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
853                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
854                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
855                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
856                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
857                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
858               END DO
859            END DO
860            !
861         ELSEIF ( ln_dynvor_een ) THEN                   !==  energy and enstrophy conserving scheme  ==!
862            DO jj = 2, jpjm1
863               DO ji = fs_2, fs_jpim1   ! vector opt.
864                  zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
865                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
866                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
867                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
868                  zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
869                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
870                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
871                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
872               END DO
873            END DO
874            !
875         ENDIF
876         !
877         ! Add tidal astronomical forcing if defined
878         IF ( ln_tide .AND. ln_tide_pot ) THEN
879            DO jj = 2, jpjm1
880               DO ji = fs_2, fs_jpim1   ! vector opt.
881                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
882                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
883                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
884                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
885               END DO
886            END DO
887         ENDIF
888         !
889         ! Add bottom stresses:
890         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:)
891         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
892         !
893         ! Add top stresses:
894         zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:)
895         zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
896         !
897         ! Surface pressure trend:
898
899         IF( ln_wd ) THEN
900           DO jj = 2, jpjm1
901              DO ji = 2, jpim1 
902                 ! Add surface pressure gradient
903                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
904                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
905                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
906                 zwy(ji,jj) = zv_spg * zcpy(ji,jj)
907              END DO
908           END DO
909         ELSE
910           DO jj = 2, jpjm1
911              DO ji = fs_2, fs_jpim1   ! vector opt.
912                 ! Add surface pressure gradient
913                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
914                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
915                 zwx(ji,jj) = zu_spg
916                 zwy(ji,jj) = zv_spg
917              END DO
918           END DO
919         END IF
920
921         !
922         ! Set next velocities:
923         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   !* Vector form
924            DO jj = 2, jpjm1
925               DO ji = fs_2, fs_jpim1   ! vector opt.
926                  ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
927                            &     + rdtbt * (                      zwx(ji,jj)   &
928                            &                                 + zu_trd(ji,jj)   &
929                            &                                 + zu_frc(ji,jj) ) & 
930                            &   ) * ssumask(ji,jj)
931
932                  va_e(ji,jj) = (                                 vn_e(ji,jj)   &
933                            &     + rdtbt * (                      zwy(ji,jj)   &
934                            &                                 + zv_trd(ji,jj)   &
935                            &                                 + zv_frc(ji,jj) ) &
936                            &   ) * ssvmask(ji,jj)
937               END DO
938            END DO
939            !
940         ELSE                                      !* Flux form
941            DO jj = 2, jpjm1
942               DO ji = fs_2, fs_jpim1   ! vector opt.
943
944                  IF( ln_wd ) THEN
945                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)
946                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)
947                  ELSE
948                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
949                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
950                  END IF
951                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
952                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
953
954                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
955                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
956                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
957                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
958                            &   ) * zhura
959
960                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
961                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
962                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
963                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
964                            &   ) * zhvra
965               END DO
966            END DO
967         ENDIF
968         !
969         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only)
970            IF( ln_wd ) THEN
971              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
972              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
973            ELSE
974              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
975              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
976            END IF
977            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
978            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
979            !
980         ENDIF
981         !                                             !* domain lateral boundary
982         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
983         !
984         !                                                 ! open boundaries
985         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
986#if defined key_agrif                                                           
987         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
988#endif
989         !                                             !* Swap
990         !                                             !  ----
991         ubb_e  (:,:) = ub_e  (:,:)
992         ub_e   (:,:) = un_e  (:,:)
993         un_e   (:,:) = ua_e  (:,:)
994         !
995         vbb_e  (:,:) = vb_e  (:,:)
996         vb_e   (:,:) = vn_e  (:,:)
997         vn_e   (:,:) = va_e  (:,:)
998         !
999         sshbb_e(:,:) = sshb_e(:,:)
1000         sshb_e (:,:) = sshn_e(:,:)
1001         sshn_e (:,:) = ssha_e(:,:)
1002
1003         !                                             !* Sum over whole bt loop
1004         !                                             !  ----------------------
1005         za1 = wgtbtp1(jn)                                   
1006         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
1007            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
1008            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
1009         ELSE                                              ! Sum transports
1010            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
1011            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
1012         ENDIF
1013         !                                   ! Sum sea level
1014         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1015         !                                                 ! ==================== !
1016      END DO                                               !        end loop      !
1017      !                                                    ! ==================== !
1018      ! -----------------------------------------------------------------------------
1019      ! Phase 3. update the general trend with the barotropic trend
1020      ! -----------------------------------------------------------------------------
1021      !
1022      ! Set advection velocity correction:
1023      zwx(:,:) = un_adv(:,:)
1024      zwy(:,:) = vn_adv(:,:)
1025      IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN     
1026         un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
1027         vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
1028      ELSE
1029         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
1030         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
1031      END IF
1032
1033      IF( ln_bt_fw ) THEN ! Save integrated transport for next computation
1034         ub2_b(:,:) = zwx(:,:)
1035         vb2_b(:,:) = zwy(:,:)
1036      ENDIF
1037      !
1038      ! Update barotropic trend:
1039      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1040         DO jk=1,jpkm1
1041            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
1042            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
1043         END DO
1044      ELSE
1045         ! At this stage, ssha has been corrected: compute new depths at velocity points
1046         DO jj = 1, jpjm1
1047            DO ji = 1, jpim1      ! NO Vector Opt.
1048               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
1049                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
1050                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1051               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
1052                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
1053                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1054            END DO
1055         END DO
1056         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1057         !
1058         DO jk=1,jpkm1
1059            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
1060            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
1061         END DO
1062         ! Save barotropic velocities not transport:
1063         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1064         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1065      ENDIF
1066      !
1067      DO jk = 1, jpkm1
1068         ! Correct velocities:
1069         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
1070         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1071         !
1072      END DO
1073      !
1074      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
1075      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic i-current
1076      !
1077#if defined key_agrif
1078      ! Save time integrated fluxes during child grid integration
1079      ! (used to update coarse grid transports at next time step)
1080      !
1081      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1082         IF( Agrif_NbStepint() == 0 ) THEN
1083            ub2_i_b(:,:) = 0._wp
1084            vb2_i_b(:,:) = 0._wp
1085         END IF
1086         !
1087         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1088         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1089         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1090      ENDIF
1091#endif     
1092      !                                   !* write time-spliting arrays in the restart
1093      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
1094      !
1095      CALL wrk_dealloc( jpi,jpj,   zsshp2_e, zhdiv )
1096      CALL wrk_dealloc( jpi,jpj,   zu_trd, zv_trd )
1097      CALL wrk_dealloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc )
1098      CALL wrk_dealloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e )
1099      CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   )
1100      CALL wrk_dealloc( jpi,jpj,   zhf )
1101      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy )
1102      !
1103      IF ( ln_diatmb ) THEN
1104         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity
1105         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity
1106      ENDIF
1107      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
1108      !
1109   END SUBROUTINE dyn_spg_ts
1110
1111
1112   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1113      !!---------------------------------------------------------------------
1114      !!                   ***  ROUTINE ts_wgt  ***
1115      !!
1116      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1117      !!----------------------------------------------------------------------
1118      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
1119      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
1120      INTEGER, INTENT(inout) :: jpit      ! cycle length   
1121      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
1122                                                         zwgt2    ! Secondary weights
1123     
1124      INTEGER ::  jic, jn, ji                      ! temporary integers
1125      REAL(wp) :: za1, za2
1126      !!----------------------------------------------------------------------
1127
1128      zwgt1(:) = 0._wp
1129      zwgt2(:) = 0._wp
1130
1131      ! Set time index when averaged value is requested
1132      IF (ll_fw) THEN
1133         jic = nn_baro
1134      ELSE
1135         jic = 2 * nn_baro
1136      ENDIF
1137
1138      ! Set primary weights:
1139      IF (ll_av) THEN
1140           ! Define simple boxcar window for primary weights
1141           ! (width = nn_baro, centered around jic)     
1142         SELECT CASE ( nn_bt_flt )
1143              CASE( 0 )  ! No averaging
1144                 zwgt1(jic) = 1._wp
1145                 jpit = jic
1146
1147              CASE( 1 )  ! Boxcar, width = nn_baro
1148                 DO jn = 1, 3*nn_baro
1149                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1150                    IF (za1 < 0.5_wp) THEN
1151                      zwgt1(jn) = 1._wp
1152                      jpit = jn
1153                    ENDIF
1154                 ENDDO
1155
1156              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
1157                 DO jn = 1, 3*nn_baro
1158                    za1 = ABS(float(jn-jic))/float(nn_baro) 
1159                    IF (za1 < 1._wp) THEN
1160                      zwgt1(jn) = 1._wp
1161                      jpit = jn
1162                    ENDIF
1163                 ENDDO
1164              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1165         END SELECT
1166
1167      ELSE ! No time averaging
1168         zwgt1(jic) = 1._wp
1169         jpit = jic
1170      ENDIF
1171   
1172      ! Set secondary weights
1173      DO jn = 1, jpit
1174        DO ji = jn, jpit
1175             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1176        END DO
1177      END DO
1178
1179      ! Normalize weigths:
1180      za1 = 1._wp / SUM(zwgt1(1:jpit))
1181      za2 = 1._wp / SUM(zwgt2(1:jpit))
1182      DO jn = 1, jpit
1183        zwgt1(jn) = zwgt1(jn) * za1
1184        zwgt2(jn) = zwgt2(jn) * za2
1185      END DO
1186      !
1187   END SUBROUTINE ts_wgt
1188
1189
1190   SUBROUTINE ts_rst( kt, cdrw )
1191      !!---------------------------------------------------------------------
1192      !!                   ***  ROUTINE ts_rst  ***
1193      !!
1194      !! ** Purpose : Read or write time-splitting arrays in restart file
1195      !!----------------------------------------------------------------------
1196      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1197      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1198      !
1199      !!----------------------------------------------------------------------
1200      !
1201      IF( TRIM(cdrw) == 'READ' ) THEN
1202         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1203         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1204         IF( .NOT.ln_bt_av ) THEN
1205            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1206            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1207            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1208            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1209            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1210            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1211         ENDIF
1212#if defined key_agrif
1213         ! Read time integrated fluxes
1214         IF ( .NOT.Agrif_Root() ) THEN
1215            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
1216            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
1217         ENDIF
1218#endif
1219      !
1220      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1221         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1222         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1223         !
1224         IF (.NOT.ln_bt_av) THEN
1225            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1226            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1227            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1228            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1229            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1230            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1231         ENDIF
1232#if defined key_agrif
1233         ! Save time integrated fluxes
1234         IF ( .NOT.Agrif_Root() ) THEN
1235            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
1236            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
1237         ENDIF
1238#endif
1239      ENDIF
1240      !
1241   END SUBROUTINE ts_rst
1242
1243
1244   SUBROUTINE dyn_spg_ts_init
1245      !!---------------------------------------------------------------------
1246      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1247      !!
1248      !! ** Purpose : Set time splitting options
1249      !!----------------------------------------------------------------------
1250      INTEGER  ::   ji ,jj              ! dummy loop indices
1251      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
1252      REAL(wp), POINTER, DIMENSION(:,:) ::   zcu
1253      !!----------------------------------------------------------------------
1254      !
1255      ! Max courant number for ext. grav. waves
1256      !
1257      CALL wrk_alloc( jpi,jpj,   zcu )
1258      !
1259      DO jj = 1, jpj
1260         DO ji =1, jpi
1261            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1262            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1263            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1264         END DO
1265      END DO
1266      !
1267      zcmax = MAXVAL( zcu(:,:) )
1268      IF( lk_mpp )   CALL mpp_max( zcmax )
1269
1270      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1271      IF( ln_bt_auto )   nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1272     
1273      rdtbt = rdt / REAL( nn_baro , wp )
1274      zcmax = zcmax * rdtbt
1275                     ! Print results
1276      IF(lwp) WRITE(numout,*)
1277      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1278      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1279      IF( ln_bt_auto ) THEN
1280         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.true. Automatically set nn_baro '
1281         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1282      ELSE
1283         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist '
1284      ENDIF
1285
1286      IF(ln_bt_av) THEN
1287         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1288      ELSE
1289         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
1290      ENDIF
1291      !
1292      !
1293      IF(ln_bt_fw) THEN
1294         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1295      ELSE
1296         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1297      ENDIF
1298      !
1299#if defined key_agrif
1300      ! Restrict the use of Agrif to the forward case only
1301      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1302#endif
1303      !
1304      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1305      SELECT CASE ( nn_bt_flt )
1306         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
1307         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1308         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1309         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1310      END SELECT
1311      !
1312      IF(lwp) WRITE(numout,*) ' '
1313      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1314      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1315      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1316      !
1317      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1318         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1319      ENDIF
1320      IF( zcmax>0.9_wp ) THEN
1321         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1322      ENDIF
1323      !
1324      CALL wrk_dealloc( jpi,jpj,   zcu )
1325      !
1326   END SUBROUTINE dyn_spg_ts_init
1327
1328   !!======================================================================
1329END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.