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

source: branches/UKMO/2015_V36_STABLE_CO6_CO5_zenv_pomsdwl/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5650

Last change on this file since 5650 was 5650, checked in by deazer, 9 years ago

Removal of svn keywords.

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 dynadv, ONLY: ln_dynadv_vec
44#if defined key_agrif
45   USE agrif_opa_interp ! agrif
46#endif
47#if defined key_asminc   
48   USE asminc          ! Assimilation increment
49#endif
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90
55   PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
56   PUBLIC dyn_spg_ts_init   !    "      "     "    "
57   PUBLIC ts_rst            !    "      "     "    "
58
59   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
60   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
61
62   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: &
63                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables
64                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables
65
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme)
69
70   ! Arrays below are saved to allow testing of the "no time averaging" option
71   ! If this option is not retained, these could be replaced by temporary arrays
72   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
73                                                   ubb_e, ub_e,     &
74                                                   vbb_e, vb_e
75
76   !! * Substitutions
77#  include "domzgr_substitute.h90"
78#  include "vectopt_loop_substitute.h90"
79   !!----------------------------------------------------------------------
80   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
81   !! $Id$
82   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
83   !!----------------------------------------------------------------------
84CONTAINS
85
86   INTEGER FUNCTION dyn_spg_ts_alloc()
87      !!----------------------------------------------------------------------
88      !!                  ***  routine dyn_spg_ts_alloc  ***
89      !!----------------------------------------------------------------------
90      INTEGER :: ierr(3)
91      !!----------------------------------------------------------------------
92      ierr(:) = 0
93
94      ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
95         &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , &
96         &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , STAT= ierr(1) )
97
98      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
99
100      IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
101                                                    &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
102
103      dyn_spg_ts_alloc = MAXVAL(ierr(:))
104
105      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
106      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
107      !
108   END FUNCTION dyn_spg_ts_alloc
109
110   SUBROUTINE dyn_spg_ts( kt )
111      !!----------------------------------------------------------------------
112      !!
113      !! ** Purpose :   
114      !!      -Compute the now trend due to the explicit time stepping
115      !!      of the quasi-linear barotropic system.
116      !!
117      !! ** Method  : 
118      !!      Barotropic variables are advanced from internal time steps
119      !!      "n"   to "n+1" if ln_bt_fw=T
120      !!      or from
121      !!      "n-1" to "n+1" if ln_bt_fw=F
122      !!      thanks to a generalized forward-backward time stepping (see ref. below).
123      !!
124      !! ** Action :
125      !!      -Update the filtered free surface at step "n+1"      : ssha
126      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
127      !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
128      !!      These are used to advect tracers and are compliant with discrete
129      !!      continuity equation taken at the baroclinic time steps. This
130      !!      ensures tracers conservation.
131      !!      -Update 3d trend (ua, va) with barotropic component.
132      !!
133      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:
134      !!              The regional oceanic modeling system (ROMS):
135      !!              a split-explicit, free-surface,
136      !!              topography-following-coordinate oceanic model.
137      !!              Ocean Modelling, 9, 347-404.
138      !!---------------------------------------------------------------------
139      !
140      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
141      !
142      LOGICAL  ::   ll_fw_start        ! if true, forward integration
143      LOGICAL  ::   ll_init             ! if true, special startup of 2d equations
144      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
145      INTEGER  ::   ikbu, ikbv, noffset      ! local integers
146      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
147      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      -
148      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      -
149      REAL(wp) ::   zu_spg, zv_spg             !   -      -
150      REAL(wp) ::   zhura, zhvra               !   -      -
151      REAL(wp) ::   za0, za1, za2, za3           !   -      -
152      !
153      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
154      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
155      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
156      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
157      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
158      REAL(wp), POINTER, DIMENSION(:,:) :: zhf
159      !!----------------------------------------------------------------------
160      !
161      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
162      !
163      !                                         !* Allocate temporary arrays
164      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
165      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  )
166      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
167      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
168      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
169      CALL wrk_alloc( jpi, jpj, zhf )
170      !
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=-2*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(:,:) + rdivisf * fwfisf(:,:) )
457      ELSE
458         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
459                &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       )
460      ENDIF
461#if defined key_asminc
462      !                                         ! Include the IAU weighted SSH increment
463      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
464         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
465      ENDIF
466#endif
467      !                                   !* Fill boundary data arrays with AGRIF
468      !                                   ! -------------------------------------
469#if defined key_agrif
470         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
471#endif
472      !
473      ! -----------------------------------------------------------------------
474      !  Phase 2 : Integration of the barotropic equations
475      ! -----------------------------------------------------------------------
476      !
477      !                                             ! ==================== !
478      !                                             !    Initialisations   !
479      !                                             ! ==================== ! 
480      ! Initialize barotropic variables:     
481      IF( ll_init )THEN
482         sshbb_e(:,:) = 0._wp
483         ubb_e  (:,:) = 0._wp
484         vbb_e  (:,:) = 0._wp
485         sshb_e (:,:) = 0._wp
486         ub_e   (:,:) = 0._wp
487         vb_e   (:,:) = 0._wp
488      ENDIF
489      !
490      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
491         sshn_e(:,:) = sshn (:,:)           
492         zun_e (:,:) = un_b (:,:)           
493         zvn_e (:,:) = vn_b (:,:)
494         !
495         hu_e  (:,:) = hu   (:,:)       
496         hv_e  (:,:) = hv   (:,:) 
497         hur_e (:,:) = hur  (:,:)   
498         hvr_e (:,:) = hvr  (:,:)
499      ELSE                                ! CENTRED integration: start from BEFORE fields
500         sshn_e(:,:) = sshb (:,:)
501         zun_e (:,:) = ub_b (:,:)         
502         zvn_e (:,:) = vb_b (:,:)
503         !
504         hu_e  (:,:) = hu_b (:,:)       
505         hv_e  (:,:) = hv_b (:,:) 
506         hur_e (:,:) = hur_b(:,:)   
507         hvr_e (:,:) = hvr_b(:,:)
508      ENDIF
509      !
510      !
511      !
512      ! Initialize sums:
513      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
514      va_b  (:,:) = 0._wp
515      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
516      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
517      zv_sum(:,:) = 0._wp
518      !                                             ! ==================== !
519      DO jn = 1, icycle                             !  sub-time-step loop  !
520         !                                          ! ==================== !
521         !                                                !* Update the forcing (BDY and tides)
522         !                                                !  ------------------
523         ! Update only tidal forcing at open boundaries
524#if defined key_tide
525         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
526         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
527#endif
528         !
529         ! Set extrapolation coefficients for predictor step:
530         IF ((jn<3).AND.ll_init) THEN      ! Forward           
531           za1 = 1._wp                                         
532           za2 = 0._wp                       
533           za3 = 0._wp                       
534         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
535           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
536           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
537           za3 =  0.281105_wp              ! za3 = bet
538         ENDIF
539
540         ! Extrapolate barotropic velocities at step jit+0.5:
541         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
542         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
543
544         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
545            !                                             !  ------------------
546            ! Extrapolate Sea Level at step jit+0.5:
547            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
548            !
549            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
550               DO ji = 2, fs_jpim1   ! Vector opt.
551                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)     &
552                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
553                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
554                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)     &
555                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
556                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
557               END DO
558            END DO
559            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
560            !
561            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
562            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
563         ELSE
564            zhup2_e (:,:) = hu(:,:)
565            zhvp2_e (:,:) = hv(:,:)
566         ENDIF
567         !                                                !* after ssh
568         !                                                !  -----------
569         ! One should enforce volume conservation at open boundaries here
570         ! considering fluxes below:
571         !
572         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
573         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
574         !
575#if defined key_agrif
576         ! Set fluxes during predictor step to ensure
577         ! volume conservation
578         IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
579            IF((nbondi == -1).OR.(nbondi == 2)) THEN
580               DO jj=1,jpj
581                  zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
582               END DO
583            ENDIF
584            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
585               DO jj=1,jpj
586                  zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
587               END DO
588            ENDIF
589            IF((nbondj == -1).OR.(nbondj == 2)) THEN
590               DO ji=1,jpi
591                  zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
592               END DO
593            ENDIF
594            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
595               DO ji=1,jpi
596                  zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
597               END DO
598            ENDIF
599         ENDIF
600#endif
601         !
602         ! Sum over sub-time-steps to compute advective velocities
603         za2 = wgtbtp2(jn)
604         zu_sum  (:,:) = zu_sum  (:,:) + za2 * zwx  (:,:) / e2u  (:,:)
605         zv_sum  (:,:) = zv_sum  (:,:) + za2 * zwy  (:,:) / e1v  (:,:)
606         !
607         ! Set next sea level:
608         DO jj = 2, jpjm1                                 
609            DO ji = fs_2, fs_jpim1   ! vector opt.
610               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
611                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e12t(ji,jj)
612            END DO
613         END DO
614         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
615         CALL lbc_lnk( ssha_e, 'T',  1._wp )
616
617#if defined key_bdy
618         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
619         IF (lk_bdy) CALL bdy_ssh( ssha_e )
620#endif
621#if defined key_agrif
622         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
623#endif
624         
625         ! Sea Surface Height at u-,v-points (vvl case only)
626         IF ( lk_vvl ) THEN                               
627            DO jj = 2, jpjm1
628               DO ji = 2, jpim1      ! NO Vector Opt.
629                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)  &
630                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
631                     &              +   e12t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
632                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)  &
633                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
634                     &              +   e12t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
635               END DO
636            END DO
637            CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
638         ENDIF   
639         !                                 
640         ! Half-step back interpolation of SSH for surface pressure computation:
641         !----------------------------------------------------------------------
642         IF ((jn==1).AND.ll_init) THEN
643           za0=1._wp                        ! Forward-backward
644           za1=0._wp                           
645           za2=0._wp
646           za3=0._wp
647         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
648           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
649           za1=-0.1666666666666_wp          ! za1 = gam
650           za2= 0.0833333333333_wp          ! za2 = eps
651           za3= 0._wp             
652         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
653           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
654           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
655           za2=0.088_wp                     ! za2 = gam
656           za3=0.013_wp                     ! za3 = eps
657         ENDIF
658
659         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
660          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
661
662         !
663         ! Compute associated depths at U and V points:
664         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN     
665            !                                       
666            DO jj = 2, jpjm1                           
667               DO ji = 2, jpim1
668                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
669                     &      * ( e12t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
670                     &      +   e12t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
671                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e12v(ji  ,jj  )  &
672                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
673                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
674                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
675                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
676               END DO
677            END DO
678         ENDIF
679         !
680         ! Add Coriolis trend:
681         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
682         ! at each time step. We however keep them constant here for optimization.
683         ! Recall that zwx and zwy arrays hold fluxes at this stage:
684         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
685         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
686         !
687         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
688            DO jj = 2, jpjm1
689               DO ji = fs_2, fs_jpim1   ! vector opt.
690                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
691                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
692                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
693                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
694                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
695                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
696               END DO
697            END DO
698            !
699         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
700            DO jj = 2, jpjm1
701               DO ji = fs_2, fs_jpim1   ! vector opt.
702                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
703                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
704                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
705                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
706                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
707                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
708               END DO
709            END DO
710            !
711         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==!
712            DO jj = 2, jpjm1
713               DO ji = fs_2, fs_jpim1   ! vector opt.
714                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
715                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
716                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
717                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
718                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
719                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
720                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
721                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
722               END DO
723            END DO
724            !
725         ENDIF
726         !
727         ! Add tidal astronomical forcing if defined
728         IF ( lk_tide.AND.ln_tide_pot ) THEN
729            DO jj = 2, jpjm1
730               DO ji = fs_2, fs_jpim1   ! vector opt.
731                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
732                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
733                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
734                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
735               END DO
736            END DO
737         ENDIF
738         !
739         ! Add bottom stresses:
740         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
741         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
742         !
743         ! Surface pressure trend:
744         DO jj = 2, jpjm1
745            DO ji = fs_2, fs_jpim1   ! vector opt.
746               ! Add surface pressure gradient
747               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
748               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
749               zwx(ji,jj) = zu_spg
750               zwy(ji,jj) = zv_spg
751            END DO
752         END DO
753         !
754         ! Set next velocities:
755         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
756            DO jj = 2, jpjm1
757               DO ji = fs_2, fs_jpim1   ! vector opt.
758                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
759                            &     + rdtbt * (                      zwx(ji,jj)   &
760                            &                                 + zu_trd(ji,jj)   &
761                            &                                 + zu_frc(ji,jj) ) & 
762                            &   ) * umask(ji,jj,1)
763
764                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
765                            &     + rdtbt * (                      zwy(ji,jj)   &
766                            &                                 + zv_trd(ji,jj)   &
767                            &                                 + zv_frc(ji,jj) ) &
768                            &   ) * vmask(ji,jj,1)
769               END DO
770            END DO
771
772         ELSE                 ! Flux form
773            DO jj = 2, jpjm1
774               DO ji = fs_2, fs_jpim1   ! vector opt.
775
776                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
777                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
778
779                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
780                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
781                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
782                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
783                            &   ) * zhura
784
785                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
786                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
787                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
788                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
789                            &   ) * zhvra
790               END DO
791            END DO
792         ENDIF
793         !
794         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
795            !                                          !  ----------------------------------------------       
796            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
797            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
798            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
799            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
800            !
801         ENDIF
802         !                                                 !* domain lateral boundary
803         !                                                 !  -----------------------
804         !
805         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
806
807#if defined key_bdy 
808                                                           ! open boundaries
809         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
810#endif
811#if defined key_agrif                                                           
812         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
813#endif
814         !                                             !* Swap
815         !                                             !  ----
816         ubb_e  (:,:) = ub_e  (:,:)
817         ub_e   (:,:) = zun_e (:,:)
818         zun_e  (:,:) = ua_e  (:,:)
819         !
820         vbb_e  (:,:) = vb_e  (:,:)
821         vb_e   (:,:) = zvn_e (:,:)
822         zvn_e  (:,:) = va_e  (:,:)
823         !
824         sshbb_e(:,:) = sshb_e(:,:)
825         sshb_e (:,:) = sshn_e(:,:)
826         sshn_e (:,:) = ssha_e(:,:)
827
828         !                                             !* Sum over whole bt loop
829         !                                             !  ----------------------
830         za1 = wgtbtp1(jn)                                   
831         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
832            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
833            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
834         ELSE                                ! Sum transports
835            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
836            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
837         ENDIF
838         !                                   ! Sum sea level
839         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
840         !                                                 ! ==================== !
841      END DO                                               !        end loop      !
842      !                                                    ! ==================== !
843      ! -----------------------------------------------------------------------------
844      ! Phase 3. update the general trend with the barotropic trend
845      ! -----------------------------------------------------------------------------
846      !
847      ! At this stage ssha holds a time averaged value
848      !                                                ! Sea Surface Height at u-,v- and f-points
849      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
850         DO jj = 1, jpjm1
851            DO ji = 1, jpim1      ! NO Vector Opt.
852               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj) &
853                  &              * ( e12t(ji  ,jj) * ssha(ji  ,jj)    &
854                  &              +   e12t(ji+1,jj) * ssha(ji+1,jj) )
855               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj) &
856                  &              * ( e12t(ji,jj  ) * ssha(ji,jj  )    &
857                  &              +   e12t(ji,jj+1) * ssha(ji,jj+1) )
858            END DO
859         END DO
860         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
861      ENDIF
862      !
863      ! Set advection velocity correction:
864      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
865         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
866         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
867      ELSE
868         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
869         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
870      END IF
871
872      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
873         ub2_b(:,:) = zu_sum(:,:)
874         vb2_b(:,:) = zv_sum(:,:)
875      ENDIF
876      !
877      ! Update barotropic trend:
878      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
879         DO jk=1,jpkm1
880            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
881            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
882         END DO
883      ELSE
884         DO jk=1,jpkm1
885            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
886            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
887         END DO
888         ! Save barotropic velocities not transport:
889         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
890         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
891      ENDIF
892      !
893      DO jk = 1, jpkm1
894         ! Correct velocities:
895         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
896         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
897         !
898      END DO
899      !
900#if defined key_agrif
901      ! Save time integrated fluxes during child grid integration
902      ! (used to update coarse grid transports)
903      ! Useless with 2nd order momentum schemes
904      !
905      IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
906         IF ( Agrif_NbStepint().EQ.0 ) THEN
907            ub2_i_b(:,:) = 0.e0
908            vb2_i_b(:,:) = 0.e0
909         END IF
910         !
911         za1 = 1._wp / REAL(Agrif_rhot(), wp)
912         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
913         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
914      ENDIF
915      !
916      !
917#endif     
918      !
919      !                                   !* write time-spliting arrays in the restart
920      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
921      !
922      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
923      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
924      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
925      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
926      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
927      CALL wrk_dealloc( jpi, jpj, zhf )
928      !
929      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
930      !
931   END SUBROUTINE dyn_spg_ts
932
933   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
934      !!---------------------------------------------------------------------
935      !!                   ***  ROUTINE ts_wgt  ***
936      !!
937      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
938      !!----------------------------------------------------------------------
939      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
940      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
941      INTEGER, INTENT(inout) :: jpit      ! cycle length   
942      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
943                                                         zwgt2    ! Secondary weights
944     
945      INTEGER ::  jic, jn, ji                      ! temporary integers
946      REAL(wp) :: za1, za2
947      !!----------------------------------------------------------------------
948
949      zwgt1(:) = 0._wp
950      zwgt2(:) = 0._wp
951
952      ! Set time index when averaged value is requested
953      IF (ll_fw) THEN
954         jic = nn_baro
955      ELSE
956         jic = 2 * nn_baro
957      ENDIF
958
959      ! Set primary weights:
960      IF (ll_av) THEN
961           ! Define simple boxcar window for primary weights
962           ! (width = nn_baro, centered around jic)     
963         SELECT CASE ( nn_bt_flt )
964              CASE( 0 )  ! No averaging
965                 zwgt1(jic) = 1._wp
966                 jpit = jic
967
968              CASE( 1 )  ! Boxcar, width = nn_baro
969                 DO jn = 1, 3*nn_baro
970                    za1 = ABS(float(jn-jic))/float(nn_baro) 
971                    IF (za1 < 0.5_wp) THEN
972                      zwgt1(jn) = 1._wp
973                      jpit = jn
974                    ENDIF
975                 ENDDO
976
977              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
978                 DO jn = 1, 3*nn_baro
979                    za1 = ABS(float(jn-jic))/float(nn_baro) 
980                    IF (za1 < 1._wp) THEN
981                      zwgt1(jn) = 1._wp
982                      jpit = jn
983                    ENDIF
984                 ENDDO
985              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
986         END SELECT
987
988      ELSE ! No time averaging
989         zwgt1(jic) = 1._wp
990         jpit = jic
991      ENDIF
992   
993      ! Set secondary weights
994      DO jn = 1, jpit
995        DO ji = jn, jpit
996             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
997        END DO
998      END DO
999
1000      ! Normalize weigths:
1001      za1 = 1._wp / SUM(zwgt1(1:jpit))
1002      za2 = 1._wp / SUM(zwgt2(1:jpit))
1003      DO jn = 1, jpit
1004        zwgt1(jn) = zwgt1(jn) * za1
1005        zwgt2(jn) = zwgt2(jn) * za2
1006      END DO
1007      !
1008   END SUBROUTINE ts_wgt
1009
1010   SUBROUTINE ts_rst( kt, cdrw )
1011      !!---------------------------------------------------------------------
1012      !!                   ***  ROUTINE ts_rst  ***
1013      !!
1014      !! ** Purpose : Read or write time-splitting arrays in restart file
1015      !!----------------------------------------------------------------------
1016      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1017      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1018      !
1019      !!----------------------------------------------------------------------
1020      !
1021      IF( TRIM(cdrw) == 'READ' ) THEN
1022         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1023         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1024         IF( .NOT.ln_bt_av ) THEN
1025            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1026            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1027            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1028            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1029            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1030            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1031         ENDIF
1032#if defined key_agrif
1033         ! Read time integrated fluxes
1034         IF ( .NOT.Agrif_Root() ) THEN
1035            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
1036            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
1037         ENDIF
1038#endif
1039      !
1040      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1041         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1042         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1043         !
1044         IF (.NOT.ln_bt_av) THEN
1045            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1046            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1047            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1048            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1049            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1050            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1051         ENDIF
1052#if defined key_agrif
1053         ! Save time integrated fluxes
1054         IF ( .NOT.Agrif_Root() ) THEN
1055            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
1056            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
1057         ENDIF
1058#endif
1059      ENDIF
1060      !
1061   END SUBROUTINE ts_rst
1062
1063   SUBROUTINE dyn_spg_ts_init( kt )
1064      !!---------------------------------------------------------------------
1065      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1066      !!
1067      !! ** Purpose : Set time splitting options
1068      !!----------------------------------------------------------------------
1069      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1070      !
1071      INTEGER  :: ji ,jj
1072      INTEGER  ::   ios                 ! Local integer output status for namelist read
1073      REAL(wp) :: zxr2, zyr2, zcmax
1074      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1075      !!
1076      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1077      &                  nn_baro, rn_bt_cmax, nn_bt_flt
1078      !!----------------------------------------------------------------------
1079      !
1080      REWIND( numnam_ref )              ! Namelist namsplit in reference namelist : time splitting parameters
1081      READ  ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1082901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1083
1084      REWIND( numnam_cfg )              ! Namelist namsplit in configuration namelist : time splitting parameters
1085      READ  ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1086902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1087      IF(lwm) WRITE ( numond, namsplit )
1088      !
1089      !         ! Max courant number for ext. grav. waves
1090      !
1091      CALL wrk_alloc( jpi, jpj, zcu )
1092      !
1093      IF (lk_vvl) THEN
1094         DO jj = 1, jpj
1095            DO ji =1, jpi
1096               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1097               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1098               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
1099            END DO
1100         END DO
1101      ELSE
1102         DO jj = 1, jpj
1103            DO ji =1, jpi
1104               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1105               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1106               zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )
1107            END DO
1108         END DO
1109      ENDIF
1110
1111      zcmax = MAXVAL(zcu(:,:))
1112      IF( lk_mpp )   CALL mpp_max( zcmax )
1113
1114      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1115      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1116     
1117      rdtbt = rdt / FLOAT(nn_baro)
1118      zcmax = zcmax * rdtbt
1119                     ! Print results
1120      IF(lwp) WRITE(numout,*)
1121      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1122      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1123      IF( ln_bt_nn_auto ) THEN
1124         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.true. Automatically set nn_baro '
1125         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1126      ELSE
1127         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1128      ENDIF
1129
1130      IF(ln_bt_av) THEN
1131         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1132      ELSE
1133         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
1134      ENDIF
1135      !
1136      !
1137      IF(ln_bt_fw) THEN
1138         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1139      ELSE
1140         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1141      ENDIF
1142      !
1143#if defined key_agrif
1144      ! Restrict the use of Agrif to the forward case only
1145      IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1146#endif
1147      !
1148      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1149      SELECT CASE ( nn_bt_flt )
1150           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac'
1151           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1152           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1153           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1154      END SELECT
1155      !
1156      IF(lwp) WRITE(numout,*) ' '
1157      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1158      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1159      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1160      !
1161      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1162         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1163      ENDIF
1164      IF ( zcmax>0.9_wp ) THEN
1165         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1166      ENDIF
1167      !
1168      CALL wrk_dealloc( jpi, jpj, zcu )
1169      !
1170   END SUBROUTINE dyn_spg_ts_init
1171
1172#else
1173   !!---------------------------------------------------------------------------
1174   !!   Default case :   Empty module   No split explicit free surface
1175   !!---------------------------------------------------------------------------
1176CONTAINS
1177   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
1178      dyn_spg_ts_alloc = 0
1179   END FUNCTION dyn_spg_ts_alloc
1180   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
1181      INTEGER, INTENT(in) :: kt
1182      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
1183   END SUBROUTINE dyn_spg_ts
1184   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
1185      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1186      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1187      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
1188   END SUBROUTINE ts_rst 
1189   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
1190      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1191      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
1192   END SUBROUTINE dyn_spg_ts_init
1193#endif
1194   
1195   !!======================================================================
1196END MODULE dynspg_ts
1197
1198
1199
Note: See TracBrowser for help on using the repository browser.