source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5910

Last change on this file since 5910 was 5910, checked in by mathiot, 5 years ago

ISF: add compatibility with time splitting

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