source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 18 months ago

Full set of changes as in the original branch

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