source: branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5605

Last change on this file since 5605 was 5605, checked in by mathiot, 6 years ago

Marine glacier termini: initial commit

  • Property svn:keywords set to Id
File size: 55.9 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
4   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
5   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
6   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
7   !!              -   ! 2008-01  (R. Benshila)  change averaging method
8   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
10   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
11   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
12   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
13   !!---------------------------------------------------------------------
14#if defined key_dynspg_ts   ||   defined key_esopa
15   !!----------------------------------------------------------------------
16   !!   'key_dynspg_ts'         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 traicw          ! ice shelf variable (fwfisf)
26   USE dynspg_oce      ! surface pressure gradient variables
27   USE phycst          ! physical constants
28   USE dynvor          ! vorticity term
29   USE bdy_par         ! for lk_bdy
30   USE bdytides        ! open boundary condition data     
31   USE bdydyn2d        ! open boundary conditions on barotropic variables
32   USE sbctide         ! tides
33   USE updtide         ! tide potential
34   USE lib_mpp         ! distributed memory computing library
35   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
36   USE prtctl          ! Print control
37   USE in_out_manager  ! I/O manager
38   USE iom             ! IOM library
39   USE restart         ! only for lrst_oce
40   USE zdf_oce         ! Bottom friction coefts
41   USE wrk_nemo        ! Memory Allocation
42   USE timing          ! Timing   
43   USE sbcapr          ! surface boundary condition: atmospheric pressure
44   USE dynadv, ONLY: ln_dynadv_vec
45#if defined key_agrif
46   USE agrif_opa_interp ! agrif
47#endif
48#if defined key_asminc   
49   USE asminc          ! Assimilation increment
50#endif
51
52   IMPLICIT NONE
53   PRIVATE
54
55   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90
56   PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
57   PUBLIC dyn_spg_ts_init   !    "      "     "    "
58   PUBLIC ts_rst            !    "      "     "    "
59
60   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
61   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
62
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: &
64                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables
65                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables
66
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme)
70
71   ! Arrays below are saved to allow testing of the "no time averaging" option
72   ! If this option is not retained, these could be replaced by temporary arrays
73   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
74                                                   ubb_e, ub_e,     &
75                                                   vbb_e, vb_e
76
77   !! * Substitutions
78#  include "domzgr_substitute.h90"
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
82   !! $Id$
83   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
84   !!----------------------------------------------------------------------
85CONTAINS
86
87   INTEGER FUNCTION dyn_spg_ts_alloc()
88      !!----------------------------------------------------------------------
89      !!                  ***  routine dyn_spg_ts_alloc  ***
90      !!----------------------------------------------------------------------
91      INTEGER :: ierr(3)
92      !!----------------------------------------------------------------------
93      ierr(:) = 0
94
95      ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
96         &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , &
97         &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , STAT= ierr(1) )
98
99      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
100
101      IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
102                                                    &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
103
104      dyn_spg_ts_alloc = MAXVAL(ierr(:))
105
106      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
107      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
108      !
109   END FUNCTION dyn_spg_ts_alloc
110
111   SUBROUTINE dyn_spg_ts( kt )
112      !!----------------------------------------------------------------------
113      !!
114      !! ** Purpose :   
115      !!      -Compute the now trend due to the explicit time stepping
116      !!      of the quasi-linear barotropic system.
117      !!
118      !! ** Method  : 
119      !!      Barotropic variables are advanced from internal time steps
120      !!      "n"   to "n+1" if ln_bt_fw=T
121      !!      or from
122      !!      "n-1" to "n+1" if ln_bt_fw=F
123      !!      thanks to a generalized forward-backward time stepping (see ref. below).
124      !!
125      !! ** Action :
126      !!      -Update the filtered free surface at step "n+1"      : ssha
127      !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
128      !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
129      !!      These are used to advect tracers and are compliant with discrete
130      !!      continuity equation taken at the baroclinic time steps. This
131      !!      ensures tracers conservation.
132      !!      -Update 3d trend (ua, va) with barotropic component.
133      !!
134      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:
135      !!              The regional oceanic modeling system (ROMS):
136      !!              a split-explicit, free-surface,
137      !!              topography-following-coordinate oceanic model.
138      !!              Ocean Modelling, 9, 347-404.
139      !!---------------------------------------------------------------------
140      !
141      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
142      !
143      LOGICAL  ::   ll_fw_start        ! if true, forward integration
144      LOGICAL  ::   ll_init             ! if true, special startup of 2d equations
145      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
146      INTEGER  ::   ikbu, ikbv, noffset      ! 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, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
248                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
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) * umask(ji,jj,1)
393             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
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      IF (ln_bt_fw) THEN                        ! Add wind forcing
423         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
424         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
425      ELSE
426         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
427         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
428      ENDIF 
429      !
430      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
431         IF (ln_bt_fw) THEN
432            DO jj = 2, jpjm1             
433               DO ji = fs_2, fs_jpim1   ! vector opt.
434                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
435                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
436                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
437                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
438               END DO
439            END DO
440         ELSE
441            DO jj = 2, jpjm1             
442               DO ji = fs_2, fs_jpim1   ! vector opt.
443                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
444                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
445                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
446                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj)
447                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
448                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
449               END DO
450            END DO
451         ENDIF
452      ENDIF
453      !                                   !* Right-Hand-Side of the barotropic ssh equation
454      !                                   ! -----------------------------------------------
455      !                                         ! Surface net water flux and rivers
456      IF (ln_bt_fw) THEN
457         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) )
458      ELSE
459         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
460                &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       &
461                                         - qtot(:,:) - qtot_b(:,:)                         )
462      ENDIF
463#if defined key_asminc
464      !                                         ! Include the IAU weighted SSH increment
465      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
466         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
467      ENDIF
468#endif
469      !                                   !* Fill boundary data arrays with AGRIF
470      !                                   ! -------------------------------------
471#if defined key_agrif
472         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
473#endif
474      !
475      ! -----------------------------------------------------------------------
476      !  Phase 2 : Integration of the barotropic equations
477      ! -----------------------------------------------------------------------
478      !
479      !                                             ! ==================== !
480      !                                             !    Initialisations   !
481      !                                             ! ==================== ! 
482      ! Initialize barotropic variables:     
483      IF( ll_init )THEN
484         sshbb_e(:,:) = 0._wp
485         ubb_e  (:,:) = 0._wp
486         vbb_e  (:,:) = 0._wp
487         sshb_e (:,:) = 0._wp
488         ub_e   (:,:) = 0._wp
489         vb_e   (:,:) = 0._wp
490      ENDIF
491      !
492      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
493         sshn_e(:,:) = sshn (:,:)           
494         zun_e (:,:) = un_b (:,:)           
495         zvn_e (:,:) = vn_b (:,:)
496         !
497         hu_e  (:,:) = hu   (:,:)       
498         hv_e  (:,:) = hv   (:,:) 
499         hur_e (:,:) = hur  (:,:)   
500         hvr_e (:,:) = hvr  (:,:)
501      ELSE                                ! CENTRED integration: start from BEFORE fields
502         sshn_e(:,:) = sshb (:,:)
503         zun_e (:,:) = ub_b (:,:)         
504         zvn_e (:,:) = vb_b (:,:)
505         !
506         hu_e  (:,:) = hu_b (:,:)       
507         hv_e  (:,:) = hv_b (:,:) 
508         hur_e (:,:) = hur_b(:,:)   
509         hvr_e (:,:) = hvr_b(:,:)
510      ENDIF
511      !
512      !
513      !
514      ! Initialize sums:
515      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
516      va_b  (:,:) = 0._wp
517      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
518      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
519      zv_sum(:,:) = 0._wp
520      !                                             ! ==================== !
521      DO jn = 1, icycle                             !  sub-time-step loop  !
522         !                                          ! ==================== !
523         !                                                !* Update the forcing (BDY and tides)
524         !                                                !  ------------------
525         ! Update only tidal forcing at open boundaries
526#if defined key_tide
527         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
528         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
529#endif
530         !
531         ! Set extrapolation coefficients for predictor step:
532         IF ((jn<3).AND.ll_init) THEN      ! Forward           
533           za1 = 1._wp                                         
534           za2 = 0._wp                       
535           za3 = 0._wp                       
536         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
537           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
538           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
539           za3 =  0.281105_wp              ! za3 = bet
540         ENDIF
541
542         ! Extrapolate barotropic velocities at step jit+0.5:
543         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
544         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
545
546         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
547            !                                             !  ------------------
548            ! Extrapolate Sea Level at step jit+0.5:
549            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
550            !
551            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
552               DO ji = 2, fs_jpim1   ! Vector opt.
553                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)     &
554                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
555                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
556                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)     &
557                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
558                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
559               END DO
560            END DO
561            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
562            !
563            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
564            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
565         ELSE
566            zhup2_e (:,:) = hu(:,:)
567            zhvp2_e (:,:) = hv(:,:)
568         ENDIF
569         !                                                !* after ssh
570         !                                                !  -----------
571         ! One should enforce volume conservation at open boundaries here
572         ! considering fluxes below:
573         !
574         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
575         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
576         !
577#if defined key_agrif
578         ! Set fluxes during predictor step to ensure
579         ! volume conservation
580         IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
581            IF((nbondi == -1).OR.(nbondi == 2)) THEN
582               DO jj=1,jpj
583                  zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
584               END DO
585            ENDIF
586            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
587               DO jj=1,jpj
588                  zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
589               END DO
590            ENDIF
591            IF((nbondj == -1).OR.(nbondj == 2)) THEN
592               DO ji=1,jpi
593                  zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
594               END DO
595            ENDIF
596            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
597               DO ji=1,jpi
598                  zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
599               END DO
600            ENDIF
601         ENDIF
602#endif
603         !
604         ! Sum over sub-time-steps to compute advective velocities
605         za2 = wgtbtp2(jn)
606         zu_sum  (:,:) = zu_sum  (:,:) + za2 * zwx  (:,:) / e2u  (:,:)
607         zv_sum  (:,:) = zv_sum  (:,:) + za2 * zwy  (:,:) / e1v  (:,:)
608         !
609         ! Set next sea level:
610         DO jj = 2, jpjm1                                 
611            DO ji = fs_2, fs_jpim1   ! vector opt.
612               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
613                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e12t(ji,jj)
614            END DO
615         END DO
616         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
617         CALL lbc_lnk( ssha_e, 'T',  1._wp )
618
619#if defined key_bdy
620         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
621         IF (lk_bdy) CALL bdy_ssh( ssha_e )
622#endif
623#if defined key_agrif
624         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
625#endif
626         
627         ! Sea Surface Height at u-,v-points (vvl case only)
628         IF ( lk_vvl ) THEN                               
629            DO jj = 2, jpjm1
630               DO ji = 2, jpim1      ! NO Vector Opt.
631                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)  &
632                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
633                     &              +   e12t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
634                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)  &
635                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
636                     &              +   e12t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
637               END DO
638            END DO
639            CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
640         ENDIF   
641         !                                 
642         ! Half-step back interpolation of SSH for surface pressure computation:
643         !----------------------------------------------------------------------
644         IF ((jn==1).AND.ll_init) THEN
645           za0=1._wp                        ! Forward-backward
646           za1=0._wp                           
647           za2=0._wp
648           za3=0._wp
649         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
650           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
651           za1=-0.1666666666666_wp          ! za1 = gam
652           za2= 0.0833333333333_wp          ! za2 = eps
653           za3= 0._wp             
654         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
655           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
656           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
657           za2=0.088_wp                     ! za2 = gam
658           za3=0.013_wp                     ! za3 = eps
659         ENDIF
660
661         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
662          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
663
664         !
665         ! Compute associated depths at U and V points:
666         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN     
667            !                                       
668            DO jj = 2, jpjm1                           
669               DO ji = 2, jpim1
670                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
671                     &      * ( e12t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
672                     &      +   e12t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
673                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e12v(ji  ,jj  )  &
674                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
675                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
676                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
677                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
678               END DO
679            END DO
680         ENDIF
681         !
682         ! Add Coriolis trend:
683         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
684         ! at each time step. We however keep them constant here for optimization.
685         ! Recall that zwx and zwy arrays hold fluxes at this stage:
686         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
687         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
688         !
689         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
690            DO jj = 2, jpjm1
691               DO ji = fs_2, fs_jpim1   ! vector opt.
692                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
693                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
694                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
695                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
696                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
697                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
698               END DO
699            END DO
700            !
701         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
702            DO jj = 2, jpjm1
703               DO ji = fs_2, fs_jpim1   ! vector opt.
704                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
705                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
706                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
707                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
708                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
709                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
710               END DO
711            END DO
712            !
713         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==!
714            DO jj = 2, jpjm1
715               DO ji = fs_2, fs_jpim1   ! vector opt.
716                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
717                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
718                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
719                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
720                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
721                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
722                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
723                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
724               END DO
725            END DO
726            !
727         ENDIF
728         !
729         ! Add tidal astronomical forcing if defined
730         IF ( lk_tide.AND.ln_tide_pot ) THEN
731            DO jj = 2, jpjm1
732               DO ji = fs_2, fs_jpim1   ! vector opt.
733                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
734                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
735                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
736                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
737               END DO
738            END DO
739         ENDIF
740         !
741         ! Add bottom stresses:
742         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
743         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
744         !
745         ! Surface pressure trend:
746         DO jj = 2, jpjm1
747            DO ji = fs_2, fs_jpim1   ! vector opt.
748               ! Add surface pressure gradient
749               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
750               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
751               zwx(ji,jj) = zu_spg
752               zwy(ji,jj) = zv_spg
753            END DO
754         END DO
755         !
756         ! Set next velocities:
757         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
758            DO jj = 2, jpjm1
759               DO ji = fs_2, fs_jpim1   ! vector opt.
760                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
761                            &     + rdtbt * (                      zwx(ji,jj)   &
762                            &                                 + zu_trd(ji,jj)   &
763                            &                                 + zu_frc(ji,jj) ) & 
764                            &   ) * umask(ji,jj,1)
765
766                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
767                            &     + rdtbt * (                      zwy(ji,jj)   &
768                            &                                 + zv_trd(ji,jj)   &
769                            &                                 + zv_frc(ji,jj) ) &
770                            &   ) * vmask(ji,jj,1)
771               END DO
772            END DO
773
774         ELSE                 ! Flux form
775            DO jj = 2, jpjm1
776               DO ji = fs_2, fs_jpim1   ! vector opt.
777
778                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
779                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
780
781                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
782                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
783                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
784                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
785                            &   ) * zhura
786
787                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
788                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
789                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
790                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
791                            &   ) * zhvra
792               END DO
793            END DO
794         ENDIF
795         !
796         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
797            !                                          !  ----------------------------------------------       
798            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
799            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
800            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
801            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
802            !
803         ENDIF
804         !                                                 !* domain lateral boundary
805         !                                                 !  -----------------------
806         !
807         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
808
809#if defined key_bdy 
810                                                           ! open boundaries
811         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
812#endif
813#if defined key_agrif                                                           
814         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
815#endif
816         !                                             !* Swap
817         !                                             !  ----
818         ubb_e  (:,:) = ub_e  (:,:)
819         ub_e   (:,:) = zun_e (:,:)
820         zun_e  (:,:) = ua_e  (:,:)
821         !
822         vbb_e  (:,:) = vb_e  (:,:)
823         vb_e   (:,:) = zvn_e (:,:)
824         zvn_e  (:,:) = va_e  (:,:)
825         !
826         sshbb_e(:,:) = sshb_e(:,:)
827         sshb_e (:,:) = sshn_e(:,:)
828         sshn_e (:,:) = ssha_e(:,:)
829
830         !                                             !* Sum over whole bt loop
831         !                                             !  ----------------------
832         za1 = wgtbtp1(jn)                                   
833         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
834            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
835            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
836         ELSE                                ! Sum transports
837            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
838            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
839         ENDIF
840         !                                   ! Sum sea level
841         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
842         !                                                 ! ==================== !
843      END DO                                               !        end loop      !
844      !                                                    ! ==================== !
845      ! -----------------------------------------------------------------------------
846      ! Phase 3. update the general trend with the barotropic trend
847      ! -----------------------------------------------------------------------------
848      !
849      ! At this stage ssha holds a time averaged value
850      !                                                ! Sea Surface Height at u-,v- and f-points
851      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
852         DO jj = 1, jpjm1
853            DO ji = 1, jpim1      ! NO Vector Opt.
854               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj) &
855                  &              * ( e12t(ji  ,jj) * ssha(ji  ,jj)    &
856                  &              +   e12t(ji+1,jj) * ssha(ji+1,jj) )
857               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj) &
858                  &              * ( e12t(ji,jj  ) * ssha(ji,jj  )    &
859                  &              +   e12t(ji,jj+1) * ssha(ji,jj+1) )
860            END DO
861         END DO
862         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
863      ENDIF
864      !
865      ! Set advection velocity correction:
866      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
867         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
868         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
869      ELSE
870         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
871         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
872      END IF
873
874      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
875         ub2_b(:,:) = zu_sum(:,:)
876         vb2_b(:,:) = zv_sum(:,:)
877      ENDIF
878      !
879      ! Update barotropic trend:
880      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
881         DO jk=1,jpkm1
882            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
883            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
884         END DO
885      ELSE
886         DO jk=1,jpkm1
887            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
888            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
889         END DO
890         ! Save barotropic velocities not transport:
891         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
892         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
893      ENDIF
894      !
895      DO jk = 1, jpkm1
896         ! Correct velocities:
897         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
898         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
899         !
900      END DO
901      !
902#if defined key_agrif
903      ! Save time integrated fluxes during child grid integration
904      ! (used to update coarse grid transports)
905      ! Useless with 2nd order momentum schemes
906      !
907      IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
908         IF ( Agrif_NbStepint().EQ.0 ) THEN
909            ub2_i_b(:,:) = 0.e0
910            vb2_i_b(:,:) = 0.e0
911         END IF
912         !
913         za1 = 1._wp / REAL(Agrif_rhot(), wp)
914         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
915         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
916      ENDIF
917      !
918      !
919#endif     
920      !
921      !                                   !* write time-spliting arrays in the restart
922      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
923      !
924      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
925      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
926      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
927      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
928      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
929      CALL wrk_dealloc( jpi, jpj, zhf )
930      !
931      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
932      !
933   END SUBROUTINE dyn_spg_ts
934
935   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
936      !!---------------------------------------------------------------------
937      !!                   ***  ROUTINE ts_wgt  ***
938      !!
939      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
940      !!----------------------------------------------------------------------
941      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
942      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
943      INTEGER, INTENT(inout) :: jpit      ! cycle length   
944      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
945                                                         zwgt2    ! Secondary weights
946     
947      INTEGER ::  jic, jn, ji                      ! temporary integers
948      REAL(wp) :: za1, za2
949      !!----------------------------------------------------------------------
950
951      zwgt1(:) = 0._wp
952      zwgt2(:) = 0._wp
953
954      ! Set time index when averaged value is requested
955      IF (ll_fw) THEN
956         jic = nn_baro
957      ELSE
958         jic = 2 * nn_baro
959      ENDIF
960
961      ! Set primary weights:
962      IF (ll_av) THEN
963           ! Define simple boxcar window for primary weights
964           ! (width = nn_baro, centered around jic)     
965         SELECT CASE ( nn_bt_flt )
966              CASE( 0 )  ! No averaging
967                 zwgt1(jic) = 1._wp
968                 jpit = jic
969
970              CASE( 1 )  ! Boxcar, width = nn_baro
971                 DO jn = 1, 3*nn_baro
972                    za1 = ABS(float(jn-jic))/float(nn_baro) 
973                    IF (za1 < 0.5_wp) THEN
974                      zwgt1(jn) = 1._wp
975                      jpit = jn
976                    ENDIF
977                 ENDDO
978
979              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
980                 DO jn = 1, 3*nn_baro
981                    za1 = ABS(float(jn-jic))/float(nn_baro) 
982                    IF (za1 < 1._wp) THEN
983                      zwgt1(jn) = 1._wp
984                      jpit = jn
985                    ENDIF
986                 ENDDO
987              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
988         END SELECT
989
990      ELSE ! No time averaging
991         zwgt1(jic) = 1._wp
992         jpit = jic
993      ENDIF
994   
995      ! Set secondary weights
996      DO jn = 1, jpit
997        DO ji = jn, jpit
998             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
999        END DO
1000      END DO
1001
1002      ! Normalize weigths:
1003      za1 = 1._wp / SUM(zwgt1(1:jpit))
1004      za2 = 1._wp / SUM(zwgt2(1:jpit))
1005      DO jn = 1, jpit
1006        zwgt1(jn) = zwgt1(jn) * za1
1007        zwgt2(jn) = zwgt2(jn) * za2
1008      END DO
1009      !
1010   END SUBROUTINE ts_wgt
1011
1012   SUBROUTINE ts_rst( kt, cdrw )
1013      !!---------------------------------------------------------------------
1014      !!                   ***  ROUTINE ts_rst  ***
1015      !!
1016      !! ** Purpose : Read or write time-splitting arrays in restart file
1017      !!----------------------------------------------------------------------
1018      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1019      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1020      !
1021      !!----------------------------------------------------------------------
1022      !
1023      IF( TRIM(cdrw) == 'READ' ) THEN
1024         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1025         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1026         IF( .NOT.ln_bt_av ) THEN
1027            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1028            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1029            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1030            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1031            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1032            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1033         ENDIF
1034#if defined key_agrif
1035         ! Read time integrated fluxes
1036         IF ( .NOT.Agrif_Root() ) THEN
1037            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
1038            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
1039         ENDIF
1040#endif
1041      !
1042      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1043         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1044         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1045         !
1046         IF (.NOT.ln_bt_av) THEN
1047            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1048            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1049            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1050            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1051            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1052            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1053         ENDIF
1054#if defined key_agrif
1055         ! Save time integrated fluxes
1056         IF ( .NOT.Agrif_Root() ) THEN
1057            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
1058            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
1059         ENDIF
1060#endif
1061      ENDIF
1062      !
1063   END SUBROUTINE ts_rst
1064
1065   SUBROUTINE dyn_spg_ts_init( kt )
1066      !!---------------------------------------------------------------------
1067      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1068      !!
1069      !! ** Purpose : Set time splitting options
1070      !!----------------------------------------------------------------------
1071      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1072      !
1073      INTEGER  :: ji ,jj
1074      INTEGER  ::   ios                 ! Local integer output status for namelist read
1075      REAL(wp) :: zxr2, zyr2, zcmax
1076      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1077      !!
1078      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1079      &                  nn_baro, rn_bt_cmax, nn_bt_flt
1080      !!----------------------------------------------------------------------
1081      !
1082      REWIND( numnam_ref )              ! Namelist namsplit in reference namelist : time splitting parameters
1083      READ  ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1084901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1085
1086      REWIND( numnam_cfg )              ! Namelist namsplit in configuration namelist : time splitting parameters
1087      READ  ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1088902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1089      IF(lwm) WRITE ( numond, namsplit )
1090      !
1091      !         ! Max courant number for ext. grav. waves
1092      !
1093      CALL wrk_alloc( jpi, jpj, zcu )
1094      !
1095      IF (lk_vvl) THEN
1096         DO jj = 1, jpj
1097            DO ji =1, jpi
1098               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1099               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1100               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
1101            END DO
1102         END DO
1103      ELSE
1104         DO jj = 1, jpj
1105            DO ji =1, jpi
1106               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1107               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1108               zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )
1109            END DO
1110         END DO
1111      ENDIF
1112
1113      zcmax = MAXVAL(zcu(:,:))
1114      IF( lk_mpp )   CALL mpp_max( zcmax )
1115
1116      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1117      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1118     
1119      rdtbt = rdt / FLOAT(nn_baro)
1120      zcmax = zcmax * rdtbt
1121                     ! Print results
1122      IF(lwp) WRITE(numout,*)
1123      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1124      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1125      IF( ln_bt_nn_auto ) THEN
1126         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.true. Automatically set nn_baro '
1127         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1128      ELSE
1129         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1130      ENDIF
1131
1132      IF(ln_bt_av) THEN
1133         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1134      ELSE
1135         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
1136      ENDIF
1137      !
1138      !
1139      IF(ln_bt_fw) THEN
1140         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1141      ELSE
1142         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1143      ENDIF
1144      !
1145#if defined key_agrif
1146      ! Restrict the use of Agrif to the forward case only
1147      IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1148#endif
1149      !
1150      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1151      SELECT CASE ( nn_bt_flt )
1152           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac'
1153           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1154           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1155           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1156      END SELECT
1157      !
1158      IF(lwp) WRITE(numout,*) ' '
1159      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1160      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1161      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1162      !
1163      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1164         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1165      ENDIF
1166      IF ( zcmax>0.9_wp ) THEN
1167         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1168      ENDIF
1169      !
1170      CALL wrk_dealloc( jpi, jpj, zcu )
1171      !
1172   END SUBROUTINE dyn_spg_ts_init
1173
1174#else
1175   !!---------------------------------------------------------------------------
1176   !!   Default case :   Empty module   No split explicit free surface
1177   !!---------------------------------------------------------------------------
1178CONTAINS
1179   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
1180      dyn_spg_ts_alloc = 0
1181   END FUNCTION dyn_spg_ts_alloc
1182   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
1183      INTEGER, INTENT(in) :: kt
1184      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
1185   END SUBROUTINE dyn_spg_ts
1186   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
1187      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1188      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1189      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
1190   END SUBROUTINE ts_rst 
1191   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
1192      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1193      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
1194   END SUBROUTINE dyn_spg_ts_init
1195#endif
1196   
1197   !!======================================================================
1198END MODULE dynspg_ts
1199
1200
1201
Note: See TracBrowser for help on using the repository browser.