source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5066

Last change on this file since 5066 was 5066, checked in by rfurner, 6 years ago

added current state of wetting and drying code to test…note it does not work

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