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

Last change on this file since 4280 was 4280, checked in by acc, 7 years ago

Branch 2013/dev_r3858_NOC_ZTC, #863. Further changes to allow PISCES SETTE tests to compile

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