source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4624

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

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.

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