New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12354

Last change on this file since 12354 was 12354, checked in by kingr, 2 years ago

Reverted change made when merging branches at r10268. Removing this line avoids introducing spurious near surface freshening when assimilating SSH obs.

File size: 56.4 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   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
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 sbcisf          ! ice shelf variable (fwfisf)
26   USE dynspg_oce      ! surface pressure gradient variables
27   USE phycst          ! physical constants
28   USE dynvor          ! vorticity term
29   USE bdy_par         ! for lk_bdy
30   USE bdytides        ! open boundary condition data     
31   USE bdydyn2d        ! open boundary conditions on barotropic variables
32   USE sbctide         ! tides
33   USE updtide         ! tide potential
34   USE sbcwave         ! surface wave
35   !
36   USE lib_mpp         ! distributed memory computing library
37   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
38   USE prtctl          ! Print control
39   USE in_out_manager  ! I/O manager
40   USE iom             ! IOM library
41   USE restart         ! only for lrst_oce
42   USE zdf_oce         ! Bottom friction coefts
43   USE wrk_nemo        ! Memory Allocation
44   USE timing          ! Timing   
45   USE sbcapr          ! surface boundary condition: atmospheric pressure
46   USE diatmb          ! Top,middle,bottom output
47   USE dynadv, ONLY: ln_dynadv_vec
48#if defined key_agrif
49   USE agrif_opa_interp ! agrif
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$
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 .or. ln_dynvor_een_old ) 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      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
146      INTEGER  ::   ikbu, ikbv, noffset      ! local integers
147      REAL(wp) ::   zmdi
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), 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      !!----------------------------------------------------------------------
162      !
163      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
164      !
165      !                                         !* Allocate temporary arrays
166      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
167      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  )
168      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
169      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
170      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
171      CALL wrk_alloc( jpi, jpj, zhf )
172      !
173      zmdi=1.e+20                               !  missing data indicator for masking
174      !                                         !* Local constant initialization
175      z1_12 = 1._wp / 12._wp 
176      z1_8  = 0.125_wp                                   
177      z1_4  = 0.25_wp
178      z1_2  = 0.5_wp     
179      zraur = 1._wp / rau0
180      !
181      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step
182        z2dt_bf = rdt
183      ELSE
184        z2dt_bf = 2.0_wp * rdt
185      ENDIF
186      z1_2dt_b = 1.0_wp / z2dt_bf 
187      !
188      ll_init = ln_bt_av                           ! if no time averaging, then no specific restart
189      ll_fw_start = .FALSE.
190      !
191                                                       ! time offset in steps for bdy data update
192      IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ;  noffset = 0 ; ENDIF
193      !
194      IF( kt == nit000 ) THEN                !* initialisation
195         !
196         IF(lwp) WRITE(numout,*)
197         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
198         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
199         IF(lwp) WRITE(numout,*)
200         !
201         IF (neuler==0) ll_init=.TRUE.
202         !
203         IF (ln_bt_fw.OR.(neuler==0)) THEN
204           ll_fw_start=.TRUE.
205           noffset = 0
206         ELSE
207           ll_fw_start=.FALSE.
208         ENDIF
209         !
210         ! Set averaging weights and cycle length:
211         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
212         !
213         !
214      ENDIF
215      !
216      ! Set arrays to remove/compute coriolis trend.
217      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
218      ! Note that these arrays are also used during barotropic loop. These are however frozen
219      ! although they should be updated in the variable volume case. Not a big approximation.
220      ! To remove this approximation, copy lines below inside barotropic loop
221      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
222      !
223      IF ( kt == nit000 .OR. lk_vvl ) THEN
224         IF ( ln_dynvor_een_old ) THEN
225            DO jj = 1, jpjm1
226               DO ji = 1, jpim1
227                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
228                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp 
229                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
230               END DO
231            END DO
232            CALL lbc_lnk( zwz, 'F', 1._wp )
233            zwz(:,:) = ff(:,:) * zwz(:,:)
234
235            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
236            DO jj = 2, jpj
237               DO ji = fs_2, jpi   ! vector opt.
238                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
239                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
240                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
241                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
242               END DO
243            END DO
244         ELSE IF ( ln_dynvor_een ) THEN
245            DO jj = 1, jpjm1
246               DO ji = 1, jpim1
247                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     &
248                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   &
249                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
250                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
251                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
252               END DO
253            END DO
254            CALL lbc_lnk( zwz, 'F', 1._wp )
255            zwz(:,:) = ff(:,:) * zwz(:,:)
256
257            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
258            DO jj = 2, jpj
259               DO ji = fs_2, jpi   ! vector opt.
260                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
261                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
262                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
263                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
264               END DO
265            END DO
266         ELSE
267            zwz(:,:) = 0._wp
268            zhf(:,:) = 0.
269            IF ( .not. ln_sco ) THEN
270! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
271!     Set it to zero for the time being
272!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
273!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
274!              ENDIF
275!              zhf(:,:) = gdepw_0(:,:,jk+1)
276            ELSE
277               zhf(:,:) = hbatf(:,:)
278            END IF
279
280            DO jj = 1, jpjm1
281               zhf(:,jj) = zhf(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
282            END DO
283
284            DO jk = 1, jpkm1
285               DO jj = 1, jpjm1
286                  zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
287               END DO
288            END DO
289            CALL lbc_lnk( zhf, 'F', 1._wp )
290            ! JC: TBC. hf should be greater than 0
291            DO jj = 1, jpj
292               DO ji = 1, jpi
293                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
294               END DO
295            END DO
296            zwz(:,:) = ff(:,:) * zwz(:,:)
297         ENDIF
298      ENDIF
299      !
300      ! If forward start at previous time step, and centered integration,
301      ! then update averaging weights:
302      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
303         ll_fw_start=.FALSE.
304         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
305      ENDIF
306                         
307      ! -----------------------------------------------------------------------------
308      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
309      ! -----------------------------------------------------------------------------
310      !     
311      !
312      !                                   !* e3*d/dt(Ua) (Vertically integrated)
313      !                                   ! --------------------------------------------------
314      zu_frc(:,:) = 0._wp
315      zv_frc(:,:) = 0._wp
316      !
317      DO jk = 1, jpkm1
318         zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
319         zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
320      END DO
321      !
322      zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
323      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
324      !
325      !
326      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
327      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
328         DO jj = 2, jpjm1
329            DO ji = fs_2, fs_jpim1   ! vector opt.
330               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
331               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
332            END DO
333         END DO
334      END DO
335             
336      !!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 
337      !!gm  Is it correct to do so ?   I think so... 
338             
339             
340      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
341      !                                   ! --------------------------------------------------------
342      zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes
343      zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
344      !
345      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
346         DO jj = 2, jpjm1
347            DO ji = fs_2, fs_jpim1   ! vector opt.
348               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
349               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
350               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
351               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
352               ! energy conserving formulation for planetary vorticity term
353               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
354               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
355            END DO
356         END DO
357         !
358      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
359         DO jj = 2, jpjm1
360            DO ji = fs_2, fs_jpim1   ! vector opt.
361               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
362                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
363               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
364                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
365               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
366               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
367            END DO
368         END DO
369         !
370      ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN  ! enstrophy and energy conserving scheme
371         DO jj = 2, jpjm1
372            DO ji = fs_2, fs_jpim1   ! vector opt.
373               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
374                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
375                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
376                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
377               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
378                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
379                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
380                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
381            END DO
382         END DO
383         !
384      ENDIF 
385      !
386      !                                   !* Right-Hand-Side of the barotropic momentum equation
387      !                                   ! ----------------------------------------------------
388      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient
389         DO jj = 2, jpjm1 
390            DO ji = fs_2, fs_jpim1   ! vector opt.
391               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
392               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
393            END DO
394         END DO
395      ENDIF
396
397      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
398         DO ji = fs_2, fs_jpim1
399             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
400             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
401          END DO
402      END DO 
403      !
404      !                 ! Add bottom stress contribution from baroclinic velocities:     
405      IF (ln_bt_fw) THEN
406         DO jj = 2, jpjm1                         
407            DO ji = fs_2, fs_jpim1   ! vector opt.
408               ikbu = mbku(ji,jj)       
409               ikbv = mbkv(ji,jj)   
410               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
411               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
412            END DO
413         END DO
414      ELSE
415         DO jj = 2, jpjm1
416            DO ji = fs_2, fs_jpim1   ! vector opt.
417               ikbu = mbku(ji,jj)       
418               ikbv = mbkv(ji,jj)   
419               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
420               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
421            END DO
422         END DO
423      ENDIF
424      !
425      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
426      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
427      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
428      !       
429      IF (ln_bt_fw) THEN                        ! Add wind forcing
430         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
431         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
432      ELSE
433         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
434         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
435      ENDIF 
436      !
437      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
438         IF (ln_bt_fw) THEN
439            DO jj = 2, jpjm1             
440               DO ji = fs_2, fs_jpim1   ! vector opt.
441                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
442                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
443                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
444                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
445               END DO
446            END DO
447         ELSE
448            DO jj = 2, jpjm1             
449               DO ji = fs_2, fs_jpim1   ! vector opt.
450                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
451                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
452                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
453                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
454                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
455                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
456               END DO
457            END DO
458         ENDIF
459      ENDIF
460      !                                   !* Right-Hand-Side of the barotropic ssh equation
461      !                                   ! -----------------------------------------------
462      !                                         ! Surface net water flux and rivers
463      IF (ln_bt_fw) THEN
464         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
465      ELSE
466         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
467                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
468      ENDIF
469      !                                   !* Fill boundary data arrays with AGRIF
470      !                                   ! -------------------------------------
471#if defined key_agrif
472         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
473#endif
474      !
475      ! -----------------------------------------------------------------------
476      !  Phase 2 : Integration of the barotropic equations
477      ! -----------------------------------------------------------------------
478      !
479      !                                             ! ==================== !
480      !                                             !    Initialisations   !
481      !                                             ! ==================== !
482      ! Initialize barotropic variables:     
483      IF( ll_init )THEN
484         sshbb_e(:,:) = 0._wp
485         ubb_e  (:,:) = 0._wp
486         vbb_e  (:,:) = 0._wp
487         sshb_e (:,:) = 0._wp
488         ub_e   (:,:) = 0._wp
489         vb_e   (:,:) = 0._wp
490      ENDIF
491      !
492      ! Moved here to allow merging with other branch 
493      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
494         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
495      ENDIF
496      !
497      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
498         sshn_e(:,:) = sshn (:,:)           
499         zun_e (:,:) = un_b (:,:)           
500         zvn_e (:,:) = vn_b (:,:)
501         !
502         hu_e  (:,:) = hu   (:,:)       
503         hv_e  (:,:) = hv   (:,:) 
504         hur_e (:,:) = hur  (:,:)   
505         hvr_e (:,:) = hvr  (:,:)
506      ELSE                                ! CENTRED integration: start from BEFORE fields
507         sshn_e(:,:) = sshb (:,:)
508         zun_e (:,:) = ub_b (:,:)         
509         zvn_e (:,:) = vb_b (:,:)
510         !
511         hu_e  (:,:) = hu_b (:,:)       
512         hv_e  (:,:) = hv_b (:,:) 
513         hur_e (:,:) = hur_b(:,:)   
514         hvr_e (:,:) = hvr_b(:,:)
515      ENDIF
516      !
517      !
518      !
519      ! Initialize sums:
520      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
521      va_b  (:,:) = 0._wp
522      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
523      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
524      zv_sum(:,:) = 0._wp
525      !                                             ! ==================== !
526      DO jn = 1, icycle                             !  sub-time-step loop  !
527         !                                          ! ==================== !
528         !                                                !* Update the forcing (BDY and tides)
529         !                                                !  ------------------
530         ! Update only tidal forcing at open boundaries
531#if defined key_tide
532         IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
533         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )
534#endif
535         !
536         ! Set extrapolation coefficients for predictor step:
537         IF ((jn<3).AND.ll_init) THEN      ! Forward           
538           za1 = 1._wp                                         
539           za2 = 0._wp                       
540           za3 = 0._wp                       
541         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
542           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
543           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
544           za3 =  0.281105_wp              ! za3 = bet
545         ENDIF
546
547         ! Extrapolate barotropic velocities at step jit+0.5:
548         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
549         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
550
551         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
552            !                                             !  ------------------
553            ! Extrapolate Sea Level at step jit+0.5:
554            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
555            !
556            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
557               DO ji = 2, fs_jpim1   ! Vector opt.
558                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)     &
559                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
560                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
561                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)     &
562                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
563                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
564               END DO
565            END DO
566            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
567            !
568            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
569            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
570         ELSE
571            zhup2_e (:,:) = hu(:,:)
572            zhvp2_e (:,:) = hv(:,:)
573         ENDIF
574         !                                                !* after ssh
575         !                                                !  -----------
576         ! One should enforce volume conservation at open boundaries here
577         ! considering fluxes below:
578         !
579         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
580         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
581         !
582#if defined key_agrif
583         ! Set fluxes during predictor step to ensure
584         ! volume conservation
585         IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
586            IF((nbondi == -1).OR.(nbondi == 2)) THEN
587               DO jj=1,jpj
588                  zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
589               END DO
590            ENDIF
591            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
592               DO jj=1,jpj
593                  zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
594               END DO
595            ENDIF
596            IF((nbondj == -1).OR.(nbondj == 2)) THEN
597               DO ji=1,jpi
598                  zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
599               END DO
600            ENDIF
601            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
602               DO ji=1,jpi
603                  zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
604               END DO
605            ENDIF
606         ENDIF
607#endif
608         !
609         ! Sum over sub-time-steps to compute advective velocities
610         za2 = wgtbtp2(jn)
611         zu_sum  (:,:) = zu_sum  (:,:) + za2 * zwx  (:,:) / e2u  (:,:)
612         zv_sum  (:,:) = zv_sum  (:,:) + za2 * zwy  (:,:) / e1v  (:,:)
613         !
614         ! Set next sea level:
615         DO jj = 2, jpjm1                                 
616            DO ji = fs_2, fs_jpim1   ! vector opt.
617               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
618                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e12t(ji,jj)
619            END DO
620         END DO
621         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
622         CALL lbc_lnk( ssha_e, 'T',  1._wp )
623
624#if defined key_bdy
625         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
626         IF (lk_bdy) CALL bdy_ssh( ssha_e )
627#endif
628#if defined key_agrif
629         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
630#endif
631         
632         ! Sea Surface Height at u-,v-points (vvl case only)
633         IF ( lk_vvl ) THEN                               
634            DO jj = 2, jpjm1
635               DO ji = 2, jpim1      ! NO Vector Opt.
636                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)  &
637                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
638                     &              +   e12t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
639                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)  &
640                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
641                     &              +   e12t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
642               END DO
643            END DO
644            CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
645         ENDIF   
646         !                                 
647         ! Half-step back interpolation of SSH for surface pressure computation:
648         !----------------------------------------------------------------------
649         IF ((jn==1).AND.ll_init) THEN
650           za0=1._wp                        ! Forward-backward
651           za1=0._wp                           
652           za2=0._wp
653           za3=0._wp
654         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
655           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
656           za1=-0.1666666666666_wp          ! za1 = gam
657           za2= 0.0833333333333_wp          ! za2 = eps
658           za3= 0._wp             
659         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
660           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
661           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
662           za2=0.088_wp                     ! za2 = gam
663           za3=0.013_wp                     ! za3 = eps
664         ENDIF
665
666         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
667          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
668
669         !
670         ! Compute associated depths at U and V points:
671         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN     
672            !                                       
673            DO jj = 2, jpjm1                           
674               DO ji = 2, jpim1
675                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
676                     &      * ( e12t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
677                     &      +   e12t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
678                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e12v(ji  ,jj  )  &
679                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
680                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
681                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
682                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
683               END DO
684            END DO
685         ENDIF
686         !
687         ! Add Coriolis trend:
688         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
689         ! at each time step. We however keep them constant here for optimization.
690         ! Recall that zwx and zwy arrays hold fluxes at this stage:
691         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
692         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
693         !
694         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
695            DO jj = 2, jpjm1
696               DO ji = fs_2, fs_jpim1   ! vector opt.
697                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
698                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
699                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
700                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
701                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
702                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
703               END DO
704            END DO
705            !
706         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
707            DO jj = 2, jpjm1
708               DO ji = fs_2, fs_jpim1   ! vector opt.
709                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
710                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
711                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
712                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
713                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
714                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
715               END DO
716            END DO
717            !
718         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==!
719            DO jj = 2, jpjm1
720               DO ji = fs_2, fs_jpim1   ! vector opt.
721                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
722                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
723                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
724                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
725                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
726                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
727                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
728                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
729               END DO
730            END DO
731            !
732         ENDIF
733         !
734         ! Add tidal astronomical forcing if defined
735         IF ( lk_tide.AND.ln_tide_pot ) THEN
736            DO jj = 2, jpjm1
737               DO ji = fs_2, fs_jpim1   ! vector opt.
738                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
739                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
740                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
741                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
742               END DO
743            END DO
744         ENDIF
745         !
746         ! Add bottom stresses:
747         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
748         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
749         !
750         ! Surface pressure trend:
751         DO jj = 2, jpjm1
752            DO ji = fs_2, fs_jpim1   ! vector opt.
753               ! Add surface pressure gradient
754               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
755               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
756               zwx(ji,jj) = zu_spg
757               zwy(ji,jj) = zv_spg
758            END DO
759         END DO
760         !
761         ! Set next velocities:
762         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
763            DO jj = 2, jpjm1
764               DO ji = fs_2, fs_jpim1   ! vector opt.
765                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
766                            &     + rdtbt * (                      zwx(ji,jj)   &
767                            &                                 + zu_trd(ji,jj)   &
768                            &                                 + zu_frc(ji,jj) ) & 
769                            &   ) * umask(ji,jj,1)
770
771                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
772                            &     + rdtbt * (                      zwy(ji,jj)   &
773                            &                                 + zv_trd(ji,jj)   &
774                            &                                 + zv_frc(ji,jj) ) &
775                            &   ) * vmask(ji,jj,1)
776               END DO
777            END DO
778
779         ELSE                 ! Flux form
780            DO jj = 2, jpjm1
781               DO ji = fs_2, fs_jpim1   ! vector opt.
782
783                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
784                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
785
786                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
787                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
788                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
789                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
790                            &   ) * zhura
791
792                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
793                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
794                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
795                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
796                            &   ) * zhvra
797               END DO
798            END DO
799         ENDIF
800         !
801         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
802            !                                          !  ----------------------------------------------       
803            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
804            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
805            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
806            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
807            !
808         ENDIF
809         !                                                 !* domain lateral boundary
810         !                                                 !  -----------------------
811         !
812         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
813
814#if defined key_bdy 
815                                                           ! open boundaries
816         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
817#endif
818#if defined key_agrif                                                           
819         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
820#endif
821         !                                             !* Swap
822         !                                             !  ----
823         ubb_e  (:,:) = ub_e  (:,:)
824         ub_e   (:,:) = zun_e (:,:)
825         zun_e  (:,:) = ua_e  (:,:)
826         !
827         vbb_e  (:,:) = vb_e  (:,:)
828         vb_e   (:,:) = zvn_e (:,:)
829         zvn_e  (:,:) = va_e  (:,:)
830         !
831         sshbb_e(:,:) = sshb_e(:,:)
832         sshb_e (:,:) = sshn_e(:,:)
833         sshn_e (:,:) = ssha_e(:,:)
834
835         !                                             !* Sum over whole bt loop
836         !                                             !  ----------------------
837         za1 = wgtbtp1(jn)                                   
838         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
839            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
840            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
841         ELSE                                ! Sum transports
842            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
843            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
844         ENDIF
845         !                                   ! Sum sea level
846         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
847         !                                                 ! ==================== !
848      END DO                                               !        end loop      !
849      !                                                    ! ==================== !
850      ! -----------------------------------------------------------------------------
851      ! Phase 3. update the general trend with the barotropic trend
852      ! -----------------------------------------------------------------------------
853      !
854      ! At this stage ssha holds a time averaged value
855      !                                                ! Sea Surface Height at u-,v- and f-points
856      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
857         DO jj = 1, jpjm1
858            DO ji = 1, jpim1      ! NO Vector Opt.
859               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj) &
860                  &              * ( e12t(ji  ,jj) * ssha(ji  ,jj)    &
861                  &              +   e12t(ji+1,jj) * ssha(ji+1,jj) )
862               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj) &
863                  &              * ( e12t(ji,jj  ) * ssha(ji,jj  )    &
864                  &              +   e12t(ji,jj+1) * ssha(ji,jj+1) )
865            END DO
866         END DO
867         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
868      ENDIF
869      !
870      ! Set advection velocity correction:
871      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
872         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
873         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
874      ELSE
875         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
876         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
877      END IF
878
879      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
880         ub2_b(:,:) = zu_sum(:,:)
881         vb2_b(:,:) = zv_sum(:,:)
882      ENDIF
883      !
884      ! Update barotropic trend:
885      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
886         DO jk=1,jpkm1
887            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
888            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
889         END DO
890      ELSE
891         DO jk=1,jpkm1
892            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
893            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
894         END DO
895         ! Save barotropic velocities not transport:
896         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
897         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
898      ENDIF
899      !
900      DO jk = 1, jpkm1
901         ! Correct velocities:
902         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
903         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
904         !
905      END DO
906      !
907#if defined key_agrif
908      ! Save time integrated fluxes during child grid integration
909      ! (used to update coarse grid transports at next time step)
910      !
911      IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
912         IF ( Agrif_NbStepint().EQ.0 ) THEN
913            ub2_i_b(:,:) = 0.e0
914            vb2_i_b(:,:) = 0.e0
915         END IF
916         !
917         za1 = 1._wp / REAL(Agrif_rhot(), wp)
918         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
919         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
920      ENDIF
921      !
922      !
923#endif     
924      !
925      !                                   !* write time-spliting arrays in the restart
926      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
927      !
928      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
929      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
930      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
931      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
932      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
933      CALL wrk_dealloc( jpi, jpj, zhf )
934      !
935      IF ( ln_diatmb ) THEN
936         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity
937         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity
938      ENDIF
939      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
940      !
941   END SUBROUTINE dyn_spg_ts
942
943   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
944      !!---------------------------------------------------------------------
945      !!                   ***  ROUTINE ts_wgt  ***
946      !!
947      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
948      !!----------------------------------------------------------------------
949      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
950      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
951      INTEGER, INTENT(inout) :: jpit      ! cycle length   
952      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
953                                                         zwgt2    ! Secondary weights
954     
955      INTEGER ::  jic, jn, ji                      ! temporary integers
956      REAL(wp) :: za1, za2
957      !!----------------------------------------------------------------------
958
959      zwgt1(:) = 0._wp
960      zwgt2(:) = 0._wp
961
962      ! Set time index when averaged value is requested
963      IF (ll_fw) THEN
964         jic = nn_baro
965      ELSE
966         jic = 2 * nn_baro
967      ENDIF
968
969      ! Set primary weights:
970      IF (ll_av) THEN
971           ! Define simple boxcar window for primary weights
972           ! (width = nn_baro, centered around jic)     
973         SELECT CASE ( nn_bt_flt )
974              CASE( 0 )  ! No averaging
975                 zwgt1(jic) = 1._wp
976                 jpit = jic
977
978              CASE( 1 )  ! Boxcar, width = nn_baro
979                 DO jn = 1, 3*nn_baro
980                    za1 = ABS(float(jn-jic))/float(nn_baro) 
981                    IF (za1 < 0.5_wp) THEN
982                      zwgt1(jn) = 1._wp
983                      jpit = jn
984                    ENDIF
985                 ENDDO
986
987              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
988                 DO jn = 1, 3*nn_baro
989                    za1 = ABS(float(jn-jic))/float(nn_baro) 
990                    IF (za1 < 1._wp) THEN
991                      zwgt1(jn) = 1._wp
992                      jpit = jn
993                    ENDIF
994                 ENDDO
995              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
996         END SELECT
997
998      ELSE ! No time averaging
999         zwgt1(jic) = 1._wp
1000         jpit = jic
1001      ENDIF
1002   
1003      ! Set secondary weights
1004      DO jn = 1, jpit
1005        DO ji = jn, jpit
1006             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1007        END DO
1008      END DO
1009
1010      ! Normalize weigths:
1011      za1 = 1._wp / SUM(zwgt1(1:jpit))
1012      za2 = 1._wp / SUM(zwgt2(1:jpit))
1013      DO jn = 1, jpit
1014        zwgt1(jn) = zwgt1(jn) * za1
1015        zwgt2(jn) = zwgt2(jn) * za2
1016      END DO
1017      !
1018   END SUBROUTINE ts_wgt
1019
1020   SUBROUTINE ts_rst( kt, cdrw )
1021      !!---------------------------------------------------------------------
1022      !!                   ***  ROUTINE ts_rst  ***
1023      !!
1024      !! ** Purpose : Read or write time-splitting arrays in restart file
1025      !!----------------------------------------------------------------------
1026      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1027      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1028      !
1029      !!----------------------------------------------------------------------
1030      !
1031      IF( TRIM(cdrw) == 'READ' ) THEN
1032         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1033         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1034         IF( .NOT.ln_bt_av ) THEN
1035            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1036            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1037            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1038            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1039            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1040            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1041         ENDIF
1042#if defined key_agrif
1043         ! Read time integrated fluxes
1044         IF ( .NOT.Agrif_Root() ) THEN
1045            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
1046            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
1047         ENDIF
1048#endif
1049      !
1050      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1051         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1052         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1053         !
1054         IF (.NOT.ln_bt_av) THEN
1055            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1056            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1057            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1058            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1059            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1060            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1061         ENDIF
1062#if defined key_agrif
1063         ! Save time integrated fluxes
1064         IF ( .NOT.Agrif_Root() ) THEN
1065            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
1066            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
1067         ENDIF
1068#endif
1069      ENDIF
1070      !
1071   END SUBROUTINE ts_rst
1072
1073   SUBROUTINE dyn_spg_ts_init( kt )
1074      !!---------------------------------------------------------------------
1075      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1076      !!
1077      !! ** Purpose : Set time splitting options
1078      !!----------------------------------------------------------------------
1079      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1080      !
1081      INTEGER  :: ji ,jj
1082      INTEGER  ::   ios                 ! Local integer output status for namelist read
1083      REAL(wp) :: zxr2, zyr2, zcmax
1084      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1085      !!
1086      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1087      &                  nn_baro, rn_bt_cmax, nn_bt_flt
1088      !!----------------------------------------------------------------------
1089      !
1090      REWIND( numnam_ref )              ! Namelist namsplit in reference namelist : time splitting parameters
1091      READ  ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1092901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1093
1094      REWIND( numnam_cfg )              ! Namelist namsplit in configuration namelist : time splitting parameters
1095      READ  ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1096902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1097      IF(lwm) WRITE ( numond, namsplit )
1098      !
1099      !         ! Max courant number for ext. grav. waves
1100      !
1101      CALL wrk_alloc( jpi, jpj, zcu )
1102      !
1103      IF (lk_vvl) THEN
1104         DO jj = 1, jpj
1105            DO ji =1, jpi
1106               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1107               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1108               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
1109            END DO
1110         END DO
1111      ELSE
1112         DO jj = 1, jpj
1113            DO ji =1, jpi
1114               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1115               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1116               zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )
1117            END DO
1118         END DO
1119      ENDIF
1120
1121      zcmax = MAXVAL(zcu(:,:))
1122      IF( lk_mpp )   CALL mpp_max( zcmax )
1123
1124      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1125      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1126     
1127      rdtbt = rdt / FLOAT(nn_baro)
1128      zcmax = zcmax * rdtbt
1129                     ! Print results
1130      IF(lwp) WRITE(numout,*)
1131      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1132      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1133      IF( ln_bt_nn_auto ) THEN
1134         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.true. Automatically set nn_baro '
1135         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1136      ELSE
1137         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1138      ENDIF
1139
1140      IF(ln_bt_av) THEN
1141         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1142      ELSE
1143         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
1144      ENDIF
1145      !
1146      !
1147      IF(ln_bt_fw) THEN
1148         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1149      ELSE
1150         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1151      ENDIF
1152      !
1153#if defined key_agrif
1154      ! Restrict the use of Agrif to the forward case only
1155      IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1156#endif
1157      !
1158      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1159      SELECT CASE ( nn_bt_flt )
1160           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac'
1161           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1162           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1163           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1164      END SELECT
1165      !
1166      IF(lwp) WRITE(numout,*) ' '
1167      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1168      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1169      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1170      !
1171      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1172         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1173      ENDIF
1174      IF ( zcmax>0.9_wp ) THEN
1175         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1176      ENDIF
1177      !
1178      CALL wrk_dealloc( jpi, jpj, zcu )
1179      !
1180   END SUBROUTINE dyn_spg_ts_init
1181
1182#else
1183   !!---------------------------------------------------------------------------
1184   !!   Default case :   Empty module   No split explicit free surface
1185   !!---------------------------------------------------------------------------
1186CONTAINS
1187   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
1188      dyn_spg_ts_alloc = 0
1189   END FUNCTION dyn_spg_ts_alloc
1190   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
1191      INTEGER, INTENT(in) :: kt
1192      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
1193   END SUBROUTINE dyn_spg_ts
1194   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
1195      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1196      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1197      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
1198   END SUBROUTINE ts_rst 
1199   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
1200      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1201      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
1202   END SUBROUTINE dyn_spg_ts_init
1203#endif
1204   
1205   !!======================================================================
1206END MODULE dynspg_ts
1207
1208
1209
Note: See TracBrowser for help on using the repository browser.