source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5014

Last change on this file since 5014 was 5014, checked in by hliu, 6 years ago

upload the modifications for W/D based on r:4826

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