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

source: branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12359

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

Removing line which subtracts SSH inc in dynspg_ts. This avoids introducing spurious near surface freshening when assimilating SSH obs.

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