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

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9178

Last change on this file since 9178 was 9178, checked in by kingr, 6 years ago

Adding changes required to apply SSH increments with VVL and/or time-splitting.

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