source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4246

Last change on this file since 4246 was 4246, checked in by acc, 8 years ago

Branch 2013/dev_r3858_NOC_ZTC, #863. Correction to divergence calculation ordering and return to single wzv routine called twice.

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