New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4370

Last change on this file since 4370 was 4370, checked in by jchanut, 10 years ago

Time split cleaning / AMM12 restart / BDY and Agrif. See tickets #1228, #1227, #1225 and #1133

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