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

source: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7606

Last change on this file since 7606 was 7606, checked in by jcastill, 6 years ago

Fix some bugs introduced due to particularities of the branch revision; make some small improvements

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