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

Last change on this file since 4288 was 4288, checked in by cbricaud, 7 years ago

correction for timesplitting

  • Property svn:keywords set to Id
File size: 55.9 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            ! JC: Simplification needed below: define ht_0 even when volume is fixed
231            IF (lk_vvl) THEN
232               zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1) 
233            ELSE
234               zht(:,:) = 0.
235               DO jk = 1, jpkm1
236                  zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
237               END DO
238            ENDIF
239
240            DO jj = 1, jpjm1
241               DO ji = 1, jpim1
242                  zwz(ji,jj) =   ( zht(ji  ,jj+1) + zht(ji+1,jj+1) +                     &
243                        &          zht(ji  ,jj  ) + zht(ji+1,jj  )   )                   &
244                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
245                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
246                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
247               END DO
248            END DO
249            CALL lbc_lnk( zwz, 'F', 1._wp )
250            zwz(:,:) = ff(:,:) * zwz(:,:)
251
252            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
253            DO jj = 2, jpj
254               DO ji = fs_2, jpi   ! vector opt.
255                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
256                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
257                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
258                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
259               END DO
260            END DO
261         ELSE
262            zwz(:,:) = 0._wp
263            zht(:,:) = 0.
264            IF ( .not. ln_sco ) THEN
265!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
266!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
267!              ENDIF
268!              zht(:,:) = gdepw_0(:,:,jk+1)
269            ELSE
270               zht(:,:) = hbatf(:,:)
271            END IF
272
273            DO jj = 1, jpjm1
274               zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1))
275            END DO
276
277            DO jk = 1, jpkm1
278               DO jj = 1, jpjm1
279                  zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
280               END DO
281            END DO
282            CALL lbc_lnk( zht, 'F', 1._wp )
283            ! JC: TBC. hf should be greater than 0
284            DO jj = 1, jpj
285               DO ji = 1, jpi
286                  IF( zht(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array
287               END DO
288            END DO
289            zwz(:,:) = ff(:,:) * zwz(:,:)
290         ENDIF
291      ENDIF
292      !
293      ! If forward start at previous time step, and centered integration,
294      ! then update averaging weights:
295      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
296         ll_fw_start=.FALSE.
297         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
298      ENDIF
299
300      ! before inverse water column height at u- and v- points
301      IF( lk_vvl ) THEN
302         zhur_b(:,:) = 0.
303         zhvr_b(:,:) = 0.
304         DO jk = 1, jpk
305            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
306            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
307         END DO
308         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) )
309         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) )
310      ELSE
311         zhur_b(:,:) = hur(:,:)
312         zhvr_b(:,:) = hvr(:,:)
313      ENDIF
314                         
315      ! -----------------------------------------------------------------------------
316      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
317      ! -----------------------------------------------------------------------------
318      !     
319      ! Some vertical sums (at now and before time steps) below could be suppressed
320      ! if one swap barotropic arrays somewhere
321      !
322      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
323      !                                   ! --------------------------------------------------
324      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp
325      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp
326      !
327      DO jk = 1, jpkm1
328#if defined key_vectopt_loop
329         DO jj = 1, 1         !Vector opt. => forced unrolling
330            DO ji = 1, jpij
331#else
332         DO jj = 1, jpj
333            DO ji = 1, jpi
334#endif
335               !        ! now trend:                                                                   
336               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
337               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)         
338               !        ! now bt transp:                   
339               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)       
340               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
341               !  ! before bt transp:
342               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
343               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
344            END DO
345         END DO
346      END DO
347      !
348      zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
349      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
350      !
351      IF( lk_vvl ) THEN
352          ub_b(:,:) = ub_b(:,:) * zhur_b(:,:)
353          vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:)
354      ELSE
355          ub_b(:,:) = ub_b(:,:) * hur(:,:)
356          vb_b(:,:) = vb_b(:,:) * hvr(:,:)
357      ENDIF
358      !
359      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
360      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
361         DO jj = 2, jpjm1
362            DO ji = fs_2, fs_jpim1   ! vector opt.
363               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
364               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
365            END DO
366         END DO
367      END DO
368      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
369      !                                   ! --------------------------------------------------------
370      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport
371      zwy(:,:) = vn_b(:,:) * e1v(:,:)
372      !
373      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
374         DO jj = 2, jpjm1
375            DO ji = fs_2, fs_jpim1   ! vector opt.
376               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
377               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
378               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
379               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
380               ! energy conserving formulation for planetary vorticity term
381               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
382               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
383            END DO
384         END DO
385         !
386      ELSEIF ( ln_dynvor_ens ) THEN             ! enstrophy conserving scheme
387         DO jj = 2, jpjm1
388            DO ji = fs_2, fs_jpim1   ! vector opt.
389               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
390                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
391               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
392                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
393               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
394               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
395            END DO
396         END DO
397         !
398      ELSEIF ( ln_dynvor_een ) THEN             ! enstrophy and energy conserving scheme
399         DO jj = 2, jpjm1
400            DO ji = fs_2, fs_jpim1   ! vector opt.
401               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
402                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
403                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
404                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
405               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
406                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
407                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
408                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
409            END DO
410         END DO
411         !
412      ENDIF 
413      !
414      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities
415      vn_b (:,:) = vn_b(:,:) * hvr(:,:) 
416      !                                   !* Right-Hand-Side of the barotropic momentum equation
417      !                                   ! ----------------------------------------------------
418      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient
419         DO jj = 2, jpjm1 
420            DO ji = fs_2, fs_jpim1   ! vector opt.
421               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
422               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
423            END DO
424         END DO
425      ENDIF
426
427      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
428         DO ji = fs_2, fs_jpim1
429             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
430             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
431          END DO
432      END DO 
433      !
434      !                 ! Add bottom stress contribution from baroclinic velocities:     
435      IF (ln_bt_fw) THEN
436         DO jj = 2, jpjm1                         
437            DO ji = fs_2, fs_jpim1   ! vector opt.
438               ikbu = mbku(ji,jj)       
439               ikbv = mbkv(ji,jj)   
440               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
441               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
442            END DO
443         END DO
444      ELSE
445         DO jj = 2, jpjm1
446            DO ji = fs_2, fs_jpim1   ! vector opt.
447               ikbu = mbku(ji,jj)       
448               ikbv = mbkv(ji,jj)   
449               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
450               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
451            END DO
452         END DO
453      ENDIF
454      !
455      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
456      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
457      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
458      !       
459      IF (ln_bt_fw) THEN                        ! Add wind forcing
460         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
461         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
462      ELSE
463         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
464         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
465      ENDIF 
466      !
467      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
468         IF (ln_bt_fw) THEN
469            DO jj = 2, jpjm1             
470               DO ji = fs_2, fs_jpim1   ! vector opt.
471                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
472                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
473                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
474                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
475               END DO
476            END DO
477         ELSE
478            DO jj = 2, jpjm1             
479               DO ji = fs_2, fs_jpim1   ! vector opt.
480                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
481                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
482                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
483                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
484                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
485                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
486               END DO
487            END DO
488         ENDIF
489      ENDIF
490      !                                   !* Right-Hand-Side of the barotropic ssh equation
491      !                                   ! -----------------------------------------------
492      !                                         ! Surface net water flux and rivers
493      IF (ln_bt_fw) THEN
494         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )
495      ELSE
496         zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:))
497      ENDIF
498#if defined key_asminc
499      !                                         ! Include the IAU weighted SSH increment
500      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
501         zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:)
502      ENDIF
503#endif
504      !
505      ! -----------------------------------------------------------------------
506      !  Phase 2 : Integration of the barotropic equations
507      ! -----------------------------------------------------------------------
508      !
509      !                                             ! ==================== !
510      !                                             !    Initialisations   !
511      !                                             ! ==================== ! 
512      ! Initialize barotropic variables:   
513      IF (ln_bt_fw) THEN                  ! FORWARD integration:  start from NOW fields                             
514         sshn_e (:,:) = sshn (:,:)           
515         zun_e  (:,:) = un_b (:,:)           
516         zvn_e  (:,:) = vn_b (:,:)
517      ELSE                                ! CENTERED integration: start from BEFORE fields
518         sshn_e (:,:) = sshb (:,:)
519         zun_e  (:,:) = ub_b (:,:)         
520         zvn_e  (:,:) = vb_b (:,:)
521      ENDIF
522      !
523      ! Initialize depths:
524      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN
525         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
526         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
527         hur_e (:,:) = zhur_b(:,:)
528         hvr_e (:,:) = zhvr_b(:,:)
529      ELSE
530         hu_e  (:,:) = hu   (:,:)       
531         hv_e  (:,:) = hv   (:,:) 
532         hur_e (:,:) = hur  (:,:)   
533         hvr_e (:,:) = hvr  (:,:)
534      ENDIF
535      !
536      IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5:
537         zhup2_e (:,:) = hu(:,:)
538         zhvp2_e (:,:) = hv(:,:)
539      ENDIF
540      !
541      ! Initialize sums:
542      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
543      va_b  (:,:) = 0._wp
544      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
545      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
546      zv_sum(:,:) = 0._wp
547      !                                             ! ==================== !
548      DO jn = 1, icycle                             !  sub-time-step loop  !
549         !                                          ! ==================== !
550         !                                                !* Update the forcing (BDY and tides)
551         !                                                !  ------------------
552         ! Update only tidal forcing at open boundaries
553#if defined key_tide
554         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
555         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
556#endif
557         !
558         ! Set extrapolation coefficients for predictor step:
559         IF ((jn<3).AND.ll_init) THEN      ! Forward           
560           za1 = 1._wp                                         
561           za2 = 0._wp                       
562           za3 = 0._wp                       
563         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
564           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
565           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
566           za3 =  0.281105_wp              ! za3 = bet
567         ENDIF
568
569         ! Extrapolate barotropic velocities at step jit+0.5:
570         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
571         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
572
573         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
574            !                                             !  ------------------
575            ! Extrapolate Sea Level at step jit+0.5:
576            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
577            !
578            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
579               DO ji = 2, fs_jpim1   ! Vector opt.
580                  zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
581                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
582                     &              +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
583                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
584                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
585                     &              +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
586               END DO
587            END DO
588            CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp )
589            !
590            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)               ! Ocean depth at U- and V-points
591            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
592         ENDIF
593         !                                                !* after ssh
594         !                                                !  -----------
595         ! One should enforce volume conservation at open boundaries here
596         ! considering fluxes below:
597         !
598         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
599         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
600         DO jj = 2, jpjm1                                 
601            DO ji = fs_2, fs_jpim1   ! vector opt.
602               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
603                  &             + zwy(ji,jj) - zwy(ji,jj-1)   &
604                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) )
605            END DO
606         END DO
607         !
608         ! Sum over sub-time-steps to compute advective velocities
609         za2 = wgtbtp2(jn)
610         zu_sum  (:,:) = zu_sum  (:,:) + za2 * ua_e  (:,:) * zhup2_e  (:,:)
611         zv_sum  (:,:) = zv_sum  (:,:) + za2 * va_e  (:,:) * zhvp2_e  (:,:)
612         !
613         ! Set next sea level:
614         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
615         CALL lbc_lnk( ssha_e, 'T',  1._wp )
616
617#if defined key_bdy
618         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
619         IF (lk_bdy) CALL bdy_ssh( ssha_e )
620#endif
621#if defined key_agrif
622         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
623#endif
624         
625         ! Sea Surface Height at u-,v-points (vvl case only)
626         IF ( lk_vvl ) THEN                               
627            DO jj = 2, jpjm1
628               DO ji = 2, jpim1      ! NO Vector Opt.
629                  zsshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 &
630                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) &
631                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) )
632                  zsshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 &
633                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) &
634                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) )
635               END DO
636            END DO
637            CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp )
638         ENDIF   
639         !                                 
640         ! Half-step back interpolation of SSH for surface pressure computation:
641         !----------------------------------------------------------------------
642         IF ((jn==1).AND.ll_init) THEN
643           za0=1._wp                        ! Forward-backward
644           za1=0._wp                           
645           za2=0._wp
646           za3=0._wp
647         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
648           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
649           za1=-0.1666666666666_wp          ! za1 = gam
650           za2= 0.0833333333333_wp          ! za2 = eps
651           za3= 0._wp             
652         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
653           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
654           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
655           za2=0.088_wp                     ! za2 = gam
656           za3=0.013_wp                     ! za3 = eps
657         ENDIF
658
659         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
660          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
661
662         !
663         ! Compute associated depths at U and V points:
664         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN     
665            !                                       
666            DO jj = 2, jpjm1                           
667               DO ji = 2, jpim1
668                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
669                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
670                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                 
671                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
672                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
673                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
674                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
675                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
676               END DO
677            END DO
678         ENDIF
679         !
680         ! Add Coriolis trend:
681         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
682         ! at each time step. We however keep them constant here for optimization.
683         ! Recall that zwx and zwy arrays hold fluxes at this stage:
684         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
685         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
686         !
687         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
688            DO jj = 2, jpjm1
689               DO ji = fs_2, fs_jpim1   ! vector opt.
690                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
691                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
692                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
693                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
694                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
695                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
696               END DO
697            END DO
698            !
699         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
700            DO jj = 2, jpjm1
701               DO ji = fs_2, fs_jpim1   ! vector opt.
702                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
703                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
704                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
705                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
706                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
707                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
708               END DO
709            END DO
710            !
711         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
712            DO jj = 2, jpjm1
713               DO ji = fs_2, fs_jpim1   ! vector opt.
714                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
715                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
716                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
717                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
718                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
719                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
720                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
721                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
722               END DO
723            END DO
724            !
725         ENDIF
726         !
727         ! Add tidal astronomical forcing if defined
728         IF ( lk_tide.AND.ln_tide_pot ) THEN
729            DO jj = 2, jpjm1
730               DO ji = fs_2, fs_jpim1   ! vector opt.
731                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
732                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
733                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
734                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
735               END DO
736            END DO
737         ENDIF
738         !
739         ! Add bottom stresses:
740         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
741         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
742         !
743         ! Surface pressure trend:
744         DO jj = 2, jpjm1
745            DO ji = fs_2, fs_jpim1   ! vector opt.
746               ! Add surface pressure gradient
747               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
748               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
749               zwx(ji,jj) = zu_spg
750               zwy(ji,jj) = zv_spg
751            END DO
752         END DO
753         !
754         ! Set next velocities:
755         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
756            DO jj = 2, jpjm1
757               DO ji = fs_2, fs_jpim1   ! vector opt.
758                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
759                            &     + rdtbt * (                      zwx(ji,jj)   &
760                            &                                 + zu_trd(ji,jj)   &
761                            &                                 + zu_frc(ji,jj) ) & 
762                            &   ) * umask(ji,jj,1)
763
764                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
765                            &     + rdtbt * (                      zwy(ji,jj)   &
766                            &                                 + zv_trd(ji,jj)   &
767                            &                                 + zv_frc(ji,jj) ) &
768                            &   ) * vmask(ji,jj,1)
769               END DO
770            END DO
771
772         ELSE                 ! Flux form
773            DO jj = 2, jpjm1
774               DO ji = fs_2, fs_jpim1   ! vector opt.
775
776                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
777                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
778
779                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
780                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
781                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
782                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
783                            &   ) * zhura
784
785                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
786                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
787                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
788                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
789                            &   ) * zhvra
790               END DO
791            END DO
792         ENDIF
793         !
794         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
795            !                                          !  ----------------------------------------------       
796            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
797            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
798            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
799            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
800            !
801         ENDIF
802         !                                                 !* domain lateral boundary
803         !                                                 !  -----------------------
804         !
805         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries
806         CALL lbc_lnk( va_e  , 'V', -1._wp )
807
808#if defined key_bdy 
809
810         pssh => ssha_e
811         phur => hur_e
812         phvr => hvr_e
813         pua2d => ua_e
814         pva2d => va_e
815         pub2d => zun_e
816         pvb2d => zvn_e
817                                     
818         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries
819#endif
820#if defined key_agrif
821         IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( kt, jn ) ! Agrif
822#endif
823         !                                             !* Swap
824         !                                             !  ----
825         ubb_e  (:,:) = ub_e  (:,:)
826         ub_e   (:,:) = zun_e (:,:)
827         zun_e  (:,:) = ua_e  (:,:)
828         !
829         vbb_e  (:,:) = vb_e  (:,:)
830         vb_e   (:,:) = zvn_e (:,:)
831         zvn_e  (:,:) = va_e  (:,:)
832         !
833         sshbb_e(:,:) = sshb_e(:,:)
834         sshb_e (:,:) = sshn_e(:,:)
835         sshn_e (:,:) = ssha_e(:,:)
836
837         !                                             !* Sum over whole bt loop
838         !                                             !  ----------------------
839         za1 = wgtbtp1(jn)                                   
840         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
841            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
842            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
843         ELSE                                ! Sum transports
844            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
845            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
846         ENDIF
847         !                                   ! Sum sea level
848         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
849         !                                                 ! ==================== !
850      END DO                                               !        end loop      !
851      !                                                    ! ==================== !
852      ! -----------------------------------------------------------------------------
853      ! Phase 3. update the general trend with the barotropic trend
854      ! -----------------------------------------------------------------------------
855      !
856      ! At this stage ssha holds a time averaged value
857      !                                                ! Sea Surface Height at u-,v- and f-points
858      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
859         DO jj = 1, jpjm1
860            DO ji = 1, jpim1      ! NO Vector Opt.
861               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
862                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
863                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
864               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
865                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
866                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
867            END DO
868         END DO
869         CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions
870      ENDIF
871      !
872      ! Set advection velocity correction:
873      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
874         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
875         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
876      ELSE
877         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
878         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
879      END IF
880
881      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
882         ub2_b(:,:) = zu_sum(:,:)
883         vb2_b(:,:) = zv_sum(:,:)
884      ENDIF
885      !
886      ! Update barotropic trend:
887      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
888         DO jk=1,jpkm1
889            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
890            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
891         END DO
892      ELSE
893         hu_e  (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
894         hv_e  (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
895         DO jk=1,jpkm1
896            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b
897            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b
898         END DO
899         ! Save barotropic velocities not transport:
900         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
901         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
902      ENDIF
903      !
904      DO jk = 1, jpkm1
905         ! Correct velocities:
906         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
907         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
908         !
909      END DO
910      !
911      !                                   !* write time-spliting arrays in the restart
912      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
913      !
914      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
915      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
916      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
917      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
918      CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b                                     )
919      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
920      CALL wrk_dealloc( jpi, jpj, zht, zhf )
921      !
922      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
923      !
924   END SUBROUTINE dyn_spg_ts
925
926   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
927      !!---------------------------------------------------------------------
928      !!                   ***  ROUTINE ts_wgt  ***
929      !!
930      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
931      !!----------------------------------------------------------------------
932      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
933      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
934      INTEGER, INTENT(inout) :: jpit      ! cycle length   
935      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
936                                                         zwgt2    ! Secondary weights
937     
938      INTEGER ::  jic, jn, ji                      ! temporary integers
939      REAL(wp) :: za1, za2
940      !!----------------------------------------------------------------------
941
942      zwgt1(:) = 0._wp
943      zwgt2(:) = 0._wp
944
945      ! Set time index when averaged value is requested
946      IF (ll_fw) THEN
947         jic = nn_baro
948      ELSE
949         jic = 2 * nn_baro
950      ENDIF
951
952      ! Set primary weights:
953      IF (ll_av) THEN
954           ! Define simple boxcar window for primary weights
955           ! (width = nn_baro, centered around jic)     
956         SELECT CASE ( nn_bt_flt )
957              CASE( 0 )  ! No averaging
958                 zwgt1(jic) = 1._wp
959                 jpit = jic
960
961              CASE( 1 )  ! Boxcar, width = nn_baro
962                 DO jn = 1, 3*nn_baro
963                    za1 = ABS(float(jn-jic))/float(nn_baro) 
964                    IF (za1 < 0.5_wp) THEN
965                      zwgt1(jn) = 1._wp
966                      jpit = jn
967                    ENDIF
968                 ENDDO
969
970              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
971                 DO jn = 1, 3*nn_baro
972                    za1 = ABS(float(jn-jic))/float(nn_baro) 
973                    IF (za1 < 1._wp) THEN
974                      zwgt1(jn) = 1._wp
975                      jpit = jn
976                    ENDIF
977                 ENDDO
978              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
979         END SELECT
980
981      ELSE ! No time averaging
982         zwgt1(jic) = 1._wp
983         jpit = jic
984      ENDIF
985   
986      ! Set secondary weights
987      DO jn = 1, jpit
988        DO ji = jn, jpit
989             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
990        END DO
991      END DO
992
993      ! Normalize weigths:
994      za1 = 1._wp / SUM(zwgt1(1:jpit))
995      za2 = 1._wp / SUM(zwgt2(1:jpit))
996      DO jn = 1, jpit
997        zwgt1(jn) = zwgt1(jn) * za1
998        zwgt2(jn) = zwgt2(jn) * za2
999      END DO
1000      !
1001   END SUBROUTINE ts_wgt
1002
1003   SUBROUTINE ts_rst( kt, cdrw )
1004      !!---------------------------------------------------------------------
1005      !!                   ***  ROUTINE ts_rst  ***
1006      !!
1007      !! ** Purpose : Read or write time-splitting arrays in restart file
1008      !!----------------------------------------------------------------------
1009      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1010      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1011      !
1012      !!----------------------------------------------------------------------
1013      !
1014      IF( TRIM(cdrw) == 'READ' ) THEN
1015         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1016         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1017         IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN
1018            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1019            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1020            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1021            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1022            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1023            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1024         ELSE
1025            sshbb_e = sshn_b                                                ! ACC GUESS WORK
1026            ubb_e   = ub_b
1027            vbb_e   = vb_b
1028            sshb_e  = sshn_b
1029            ub_e    = ub_b
1030            vb_e    = vb_b
1031         ENDIF
1032      !
1033      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1034         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1035         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1036         !
1037         IF (.NOT.ln_bt_av) THEN
1038            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1039            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1040            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1041            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1042            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1043            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1044         ENDIF
1045      ENDIF
1046      !
1047   END SUBROUTINE ts_rst
1048
1049   SUBROUTINE dyn_spg_ts_init( kt )
1050      !!---------------------------------------------------------------------
1051      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1052      !!
1053      !! ** Purpose : Set time splitting options
1054      !!----------------------------------------------------------------------
1055      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1056      !
1057      INTEGER  :: ji ,jj, jk
1058      REAL(wp) :: zxr2, zyr2, zcmax
1059      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht
1060      !!
1061!      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1062!      &                  nn_baro, rn_bt_cmax, nn_bt_flt
1063      !!----------------------------------------------------------------------
1064!      REWIND( numnam )              !* Namelist namsplit: split-explicit free surface
1065!      READ  ( numnam, namsplit )
1066      !         ! Max courant number for ext. grav. waves
1067      !
1068      CALL wrk_alloc( jpi, jpj, zcu, zht )
1069      !
1070      ! JC: Simplification needed below: define ht_0 even when volume is fixed
1071      IF (lk_vvl) THEN
1072         zht(:,:) = ht_0(:,:) * tmask(:,:,1)
1073      ELSE
1074         zht(:,:) = 0.
1075         DO jk = 1, jpkm1
1076            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
1077         END DO
1078      ENDIF
1079
1080      DO jj = 1, jpj
1081         DO ji =1, jpi
1082            zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1083            zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1084            zcu(ji,jj) = sqrt(grav*zht(ji,jj)*(zxr2 + zyr2) )
1085         END DO
1086      END DO
1087
1088      zcmax = MAXVAL(zcu(:,:))
1089      IF( lk_mpp )   CALL mpp_max( zcmax )
1090
1091      ! Estimate number of iterations to satisfy a max courant number=0.8
1092      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1093     
1094      rdtbt = rdt / FLOAT(nn_baro)
1095      zcmax = zcmax * rdtbt
1096                     ! Print results
1097      IF(lwp) WRITE(numout,*)
1098      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1099      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1100      IF( ln_bt_nn_auto ) THEN
1101         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
1102         IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1103      ELSE
1104         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1105      ENDIF
1106      IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1107      IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1108      IF(lwp) WRITE(numout,*) ' Maximum Courant number is   :', zcmax
1109
1110      IF(ln_bt_av) THEN
1111         IF(lwp) WRITE(numout,*) ' ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1112      ELSE
1113         IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1114      ENDIF
1115      !
1116      !
1117      IF(ln_bt_fw) THEN
1118         IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables '
1119      ELSE
1120         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables '
1121      ENDIF
1122      !
1123      IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1124      SELECT CASE ( nn_bt_flt )
1125           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '      Dirac'
1126           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = nn_baro'
1127           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = 2*nn_baro' 
1128           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1129      END SELECT
1130      !
1131      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1132         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1133      ENDIF
1134      IF ( zcmax>0.9_wp ) THEN
1135         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1136      ENDIF
1137      !
1138      CALL wrk_dealloc( jpi, jpj, zcu, zht )
1139      !
1140   END SUBROUTINE dyn_spg_ts_init
1141
1142#else
1143   !!---------------------------------------------------------------------------
1144   !!   Default case :   Empty module   No standard free surface constant volume
1145   !!---------------------------------------------------------------------------
1146
1147   USE par_kind
1148   LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping
1149CONTAINS
1150   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
1151      dyn_spg_ts_alloc = 0
1152   END FUNCTION dyn_spg_ts_alloc
1153   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
1154      INTEGER, INTENT(in) :: kt
1155      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
1156   END SUBROUTINE dyn_spg_ts
1157   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
1158      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1159      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1160      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
1161   END SUBROUTINE ts_rst 
1162   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
1163      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1164      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
1165   END SUBROUTINE dyn_spg_ts_init
1166#endif
1167   
1168   !!======================================================================
1169END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.