source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 19 months ago

GMED 450 add flush after prints

File size: 56.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      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
147      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      -
148      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      -
149      REAL(wp) ::   zu_spg, zv_spg             !   -      -
150      REAL(wp) ::   zhura, zhvra               !   -      -
151      REAL(wp) ::   za0, za1, za2, za3           !   -      -
152      !
153      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
154      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
155      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv
156      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
157      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
158      REAL(wp), POINTER, DIMENSION(:,:) :: zhf
159      !!----------------------------------------------------------------------
160      !
161      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
162      !
163      !                                         !* Allocate temporary arrays
164      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv )
165      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  )
166      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc)
167      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
168      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
169      CALL wrk_alloc( jpi, jpj, zhf )
170      !
171      !                                         !* Local constant initialization
172      z1_12 = 1._wp / 12._wp 
173      z1_8  = 0.125_wp                                   
174      z1_4  = 0.25_wp
175      z1_2  = 0.5_wp     
176      zraur = 1._wp / rau0
177      !
178      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step
179        z2dt_bf = rdt
180      ELSE
181        z2dt_bf = 2.0_wp * rdt
182      ENDIF
183      z1_2dt_b = 1.0_wp / z2dt_bf 
184      !
185      ll_init = ln_bt_av                           ! if no time averaging, then no specific restart
186      ll_fw_start = .FALSE.
187      !
188                                                       ! time offset in steps for bdy data update
189      IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ;  noffset = 0 ; ENDIF
190      !
191      IF( kt == nit000 ) THEN                !* initialisation
192         !
193         IF(lwp) WRITE(numout,*)
194         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
195         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
196         IF(lwp) WRITE(numout,*)
197         IF(lflush) CALL flush(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(:,:) + fwfisf(:,:) )
458      ELSE
459         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
460                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
461      ENDIF
462#if defined key_asminc
463      !                                         ! Include the IAU weighted SSH increment
464      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
465         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
466      ENDIF
467#endif
468      !                                   !* Fill boundary data arrays for AGRIF
469      !                                   ! ------------------------------------
470#if defined key_agrif
471         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
472#endif
473      !
474      ! -----------------------------------------------------------------------
475      !  Phase 2 : Integration of the barotropic equations
476      ! -----------------------------------------------------------------------
477      !
478      !                                             ! ==================== !
479      !                                             !    Initialisations   !
480      !                                             ! ==================== ! 
481      ! Initialize barotropic variables:     
482      IF( ll_init )THEN
483         sshbb_e(:,:) = 0._wp
484         ubb_e  (:,:) = 0._wp
485         vbb_e  (:,:) = 0._wp
486         sshb_e (:,:) = 0._wp
487         ub_e   (:,:) = 0._wp
488         vb_e   (:,:) = 0._wp
489      ENDIF
490      !
491      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
492         sshn_e(:,:) = sshn (:,:)           
493         zun_e (:,:) = un_b (:,:)           
494         zvn_e (:,:) = vn_b (:,:)
495         !
496         hu_e  (:,:) = hu   (:,:)       
497         hv_e  (:,:) = hv   (:,:) 
498         hur_e (:,:) = hur  (:,:)   
499         hvr_e (:,:) = hvr  (:,:)
500      ELSE                                ! CENTRED integration: start from BEFORE fields
501         sshn_e(:,:) = sshb (:,:)
502         zun_e (:,:) = ub_b (:,:)         
503         zvn_e (:,:) = vb_b (:,:)
504         !
505         hu_e  (:,:) = hu_b (:,:)       
506         hv_e  (:,:) = hv_b (:,:) 
507         hur_e (:,:) = hur_b(:,:)   
508         hvr_e (:,:) = hvr_b(:,:)
509      ENDIF
510      !
511      !
512      !
513      ! Initialize sums:
514      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)         
515      va_b  (:,:) = 0._wp
516      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
517      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop
518      zv_sum(:,:) = 0._wp
519      !                                             ! ==================== !
520      DO jn = 1, icycle                             !  sub-time-step loop  !
521         !                                          ! ==================== !
522         !                                                !* Update the forcing (BDY and tides)
523         !                                                !  ------------------
524         ! Update only tidal forcing at open boundaries
525#if defined key_tide
526         IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
527         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )
528#endif
529         !
530         ! Set extrapolation coefficients for predictor step:
531         IF ((jn<3).AND.ll_init) THEN      ! Forward           
532           za1 = 1._wp                                         
533           za2 = 0._wp                       
534           za3 = 0._wp                       
535         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
536           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
537           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
538           za3 =  0.281105_wp              ! za3 = bet
539         ENDIF
540
541         ! Extrapolate barotropic velocities at step jit+0.5:
542         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
543         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
544
545         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
546            !                                             !  ------------------
547            ! Extrapolate Sea Level at step jit+0.5:
548            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
549            !
550            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
551               DO ji = 2, fs_jpim1   ! Vector opt.
552                  zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)     &
553                     &              * ( e12t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
554                     &              +   e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
555                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)     &
556                     &              * ( e12t(ji,jj  ) * zsshp2_e(ji,jj  )  &
557                     &              +   e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
558               END DO
559            END DO
560            CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
561            !
562            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
563            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
564         ELSE
565            zhup2_e (:,:) = hu(:,:)
566            zhvp2_e (:,:) = hv(:,:)
567         ENDIF
568         !                                                !* after ssh
569         !                                                !  -----------
570         ! One should enforce volume conservation at open boundaries here
571         ! considering fluxes below:
572         !
573         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
574         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
575         !
576#if defined key_agrif
577         ! Set fluxes during predictor step to ensure
578         ! volume conservation
579         IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
580            IF((nbondi == -1).OR.(nbondi == 2)) THEN
581               DO jj=1,jpj
582                  zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
583               END DO
584            ENDIF
585            IF((nbondi ==  1).OR.(nbondi == 2)) THEN
586               DO jj=1,jpj
587                  zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
588               END DO
589            ENDIF
590            IF((nbondj == -1).OR.(nbondj == 2)) THEN
591               DO ji=1,jpi
592                  zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
593               END DO
594            ENDIF
595            IF((nbondj ==  1).OR.(nbondj == 2)) THEN
596               DO ji=1,jpi
597                  zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
598               END DO
599            ENDIF
600         ENDIF
601#endif
602         !
603         ! Sum over sub-time-steps to compute advective velocities
604         za2 = wgtbtp2(jn)
605         zu_sum  (:,:) = zu_sum  (:,:) + za2 * zwx  (:,:) / e2u  (:,:)
606         zv_sum  (:,:) = zv_sum  (:,:) + za2 * zwy  (:,:) / e1v  (:,:)
607         !
608         ! Set next sea level:
609         DO jj = 2, jpjm1                                 
610            DO ji = fs_2, fs_jpim1   ! vector opt.
611               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
612                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e12t(ji,jj)
613            END DO
614         END DO
615         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
616         CALL lbc_lnk( ssha_e, 'T',  1._wp )
617
618#if defined key_bdy
619         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
620         IF (lk_bdy) CALL bdy_ssh( ssha_e )
621#endif
622#if defined key_agrif
623         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
624#endif
625         
626         ! Sea Surface Height at u-,v-points (vvl case only)
627         IF ( lk_vvl ) THEN                               
628            DO jj = 2, jpjm1
629               DO ji = 2, jpim1      ! NO Vector Opt.
630                  zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj)  &
631                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
632                     &              +   e12t(ji+1,jj  ) * ssha_e(ji+1,jj  ) )
633                  zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj)  &
634                     &              * ( e12t(ji  ,jj  ) * ssha_e(ji  ,jj  ) &
635                     &              +   e12t(ji  ,jj+1) * ssha_e(ji  ,jj+1) )
636               END DO
637            END DO
638            CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
639         ENDIF   
640         !                                 
641         ! Half-step back interpolation of SSH for surface pressure computation:
642         !----------------------------------------------------------------------
643         IF ((jn==1).AND.ll_init) THEN
644           za0=1._wp                        ! Forward-backward
645           za1=0._wp                           
646           za2=0._wp
647           za3=0._wp
648         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
649           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
650           za1=-0.1666666666666_wp          ! za1 = gam
651           za2= 0.0833333333333_wp          ! za2 = eps
652           za3= 0._wp             
653         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
654           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps   
655           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
656           za2=0.088_wp                     ! za2 = gam
657           za3=0.013_wp                     ! za3 = eps
658         ENDIF
659
660         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
661          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
662
663         !
664         ! Compute associated depths at U and V points:
665         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN     
666            !                                       
667            DO jj = 2, jpjm1                           
668               DO ji = 2, jpim1
669                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
670                     &      * ( e12t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
671                     &      +   e12t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
672                  zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e12v(ji  ,jj  )  &
673                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
674                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
675                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
676                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
677               END DO
678            END DO
679         ENDIF
680         !
681         ! Add Coriolis trend:
682         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated
683         ! at each time step. We however keep them constant here for optimization.
684         ! Recall that zwx and zwy arrays hold fluxes at this stage:
685         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
686         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
687         !
688         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
689            DO jj = 2, jpjm1
690               DO ji = fs_2, fs_jpim1   ! vector opt.
691                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
692                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
693                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
694                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
695                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
696                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
697               END DO
698            END DO
699            !
700         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
701            DO jj = 2, jpjm1
702               DO ji = fs_2, fs_jpim1   ! vector opt.
703                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
704                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
705                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
706                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
707                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
708                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
709               END DO
710            END DO
711            !
712         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==!
713            DO jj = 2, jpjm1
714               DO ji = fs_2, fs_jpim1   ! vector opt.
715                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
716                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
717                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
718                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
719                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
720                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
721                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
722                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
723               END DO
724            END DO
725            !
726         ENDIF
727         !
728         ! Add tidal astronomical forcing if defined
729         IF ( lk_tide.AND.ln_tide_pot ) THEN
730            DO jj = 2, jpjm1
731               DO ji = fs_2, fs_jpim1   ! vector opt.
732                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
733                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
734                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
735                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
736               END DO
737            END DO
738         ENDIF
739         !
740         ! Add bottom stresses:
741         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
742         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
743         !
744         ! Surface pressure trend:
745         DO jj = 2, jpjm1
746            DO ji = fs_2, fs_jpim1   ! vector opt.
747               ! Add surface pressure gradient
748               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
749               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
750               zwx(ji,jj) = zu_spg
751               zwy(ji,jj) = zv_spg
752            END DO
753         END DO
754         !
755         ! Set next velocities:
756         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
757            DO jj = 2, jpjm1
758               DO ji = fs_2, fs_jpim1   ! vector opt.
759                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
760                            &     + rdtbt * (                      zwx(ji,jj)   &
761                            &                                 + zu_trd(ji,jj)   &
762                            &                                 + zu_frc(ji,jj) ) & 
763                            &   ) * umask(ji,jj,1)
764
765                  va_e(ji,jj) = (                                zvn_e(ji,jj)   &
766                            &     + rdtbt * (                      zwy(ji,jj)   &
767                            &                                 + zv_trd(ji,jj)   &
768                            &                                 + zv_frc(ji,jj) ) &
769                            &   ) * vmask(ji,jj,1)
770               END DO
771            END DO
772
773         ELSE                 ! Flux form
774            DO jj = 2, jpjm1
775               DO ji = fs_2, fs_jpim1   ! vector opt.
776
777                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
778                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
779
780                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
781                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
782                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
783                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
784                            &   ) * zhura
785
786                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   &
787                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
788                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
789                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
790                            &   ) * zhvra
791               END DO
792            END DO
793         ENDIF
794         !
795         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
796            !                                          !  ----------------------------------------------       
797            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
798            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
799            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
800            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
801            !
802         ENDIF
803         !                                                 !* domain lateral boundary
804         !                                                 !  -----------------------
805         !
806         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
807
808#if defined key_bdy 
809                                                           ! open boundaries
810         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
811#endif
812#if defined key_agrif                                                           
813         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
814#endif
815         !                                             !* Swap
816         !                                             !  ----
817         ubb_e  (:,:) = ub_e  (:,:)
818         ub_e   (:,:) = zun_e (:,:)
819         zun_e  (:,:) = ua_e  (:,:)
820         !
821         vbb_e  (:,:) = vb_e  (:,:)
822         vb_e   (:,:) = zvn_e (:,:)
823         zvn_e  (:,:) = va_e  (:,:)
824         !
825         sshbb_e(:,:) = sshb_e(:,:)
826         sshb_e (:,:) = sshn_e(:,:)
827         sshn_e (:,:) = ssha_e(:,:)
828
829         !                                             !* Sum over whole bt loop
830         !                                             !  ----------------------
831         za1 = wgtbtp1(jn)                                   
832         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
833            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
834            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
835         ELSE                                ! Sum transports
836            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
837            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
838         ENDIF
839         !                                   ! Sum sea level
840         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
841         !                                                 ! ==================== !
842      END DO                                               !        end loop      !
843      !                                                    ! ==================== !
844      ! -----------------------------------------------------------------------------
845      ! Phase 3. update the general trend with the barotropic trend
846      ! -----------------------------------------------------------------------------
847      !
848      ! At this stage ssha holds a time averaged value
849      !                                                ! Sea Surface Height at u-,v- and f-points
850      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
851         DO jj = 1, jpjm1
852            DO ji = 1, jpim1      ! NO Vector Opt.
853               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e12u(ji,jj) &
854                  &              * ( e12t(ji  ,jj) * ssha(ji  ,jj)    &
855                  &              +   e12t(ji+1,jj) * ssha(ji+1,jj) )
856               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e12v(ji,jj) &
857                  &              * ( e12t(ji,jj  ) * ssha(ji,jj  )    &
858                  &              +   e12t(ji,jj+1) * ssha(ji,jj+1) )
859            END DO
860         END DO
861         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
862      ENDIF
863      !
864      ! Set advection velocity correction:
865      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
866         un_adv(:,:) = zu_sum(:,:)*hur(:,:)
867         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
868      ELSE
869         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
870         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
871      END IF
872
873      IF (ln_bt_fw) THEN ! Save integrated transport for next computation
874         ub2_b(:,:) = zu_sum(:,:)
875         vb2_b(:,:) = zv_sum(:,:)
876      ENDIF
877      !
878      ! Update barotropic trend:
879      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
880         DO jk=1,jpkm1
881            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
882            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
883         END DO
884      ELSE
885         DO jk=1,jpkm1
886            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
887            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
888         END DO
889         ! Save barotropic velocities not transport:
890         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
891         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
892      ENDIF
893      !
894      DO jk = 1, jpkm1
895         ! Correct velocities:
896         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
897         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
898         !
899      END DO
900      !
901#if defined key_agrif
902      ! Save time integrated fluxes during child grid integration
903      ! (used to update coarse grid transports at next time step)
904      !
905      IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
906         IF ( Agrif_NbStepint().EQ.0 ) THEN
907            ub2_i_b(:,:) = 0.e0
908            vb2_i_b(:,:) = 0.e0
909         END IF
910         !
911         za1 = 1._wp / REAL(Agrif_rhot(), wp)
912         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
913         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
914      ENDIF
915      !
916      !
917#endif     
918      !
919      !                                   !* write time-spliting arrays in the restart
920      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' )
921      !
922      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
923      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
924      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
925      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
926      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
927      CALL wrk_dealloc( jpi, jpj, zhf )
928      !
929      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
930      !
931   END SUBROUTINE dyn_spg_ts
932
933   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
934      !!---------------------------------------------------------------------
935      !!                   ***  ROUTINE ts_wgt  ***
936      !!
937      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
938      !!----------------------------------------------------------------------
939      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
940      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
941      INTEGER, INTENT(inout) :: jpit      ! cycle length   
942      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
943                                                         zwgt2    ! Secondary weights
944     
945      INTEGER ::  jic, jn, ji                      ! temporary integers
946      REAL(wp) :: za1, za2
947      !!----------------------------------------------------------------------
948
949      zwgt1(:) = 0._wp
950      zwgt2(:) = 0._wp
951
952      ! Set time index when averaged value is requested
953      IF (ll_fw) THEN
954         jic = nn_baro
955      ELSE
956         jic = 2 * nn_baro
957      ENDIF
958
959      ! Set primary weights:
960      IF (ll_av) THEN
961           ! Define simple boxcar window for primary weights
962           ! (width = nn_baro, centered around jic)     
963         SELECT CASE ( nn_bt_flt )
964              CASE( 0 )  ! No averaging
965                 zwgt1(jic) = 1._wp
966                 jpit = jic
967
968              CASE( 1 )  ! Boxcar, width = nn_baro
969                 DO jn = 1, 3*nn_baro
970                    za1 = ABS(float(jn-jic))/float(nn_baro) 
971                    IF (za1 < 0.5_wp) THEN
972                      zwgt1(jn) = 1._wp
973                      jpit = jn
974                    ENDIF
975                 ENDDO
976
977              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
978                 DO jn = 1, 3*nn_baro
979                    za1 = ABS(float(jn-jic))/float(nn_baro) 
980                    IF (za1 < 1._wp) THEN
981                      zwgt1(jn) = 1._wp
982                      jpit = jn
983                    ENDIF
984                 ENDDO
985              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
986         END SELECT
987
988      ELSE ! No time averaging
989         zwgt1(jic) = 1._wp
990         jpit = jic
991      ENDIF
992   
993      ! Set secondary weights
994      DO jn = 1, jpit
995        DO ji = jn, jpit
996             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
997        END DO
998      END DO
999
1000      ! Normalize weigths:
1001      za1 = 1._wp / SUM(zwgt1(1:jpit))
1002      za2 = 1._wp / SUM(zwgt2(1:jpit))
1003      DO jn = 1, jpit
1004        zwgt1(jn) = zwgt1(jn) * za1
1005        zwgt2(jn) = zwgt2(jn) * za2
1006      END DO
1007      !
1008   END SUBROUTINE ts_wgt
1009
1010   SUBROUTINE ts_rst( kt, cdrw )
1011      !!---------------------------------------------------------------------
1012      !!                   ***  ROUTINE ts_rst  ***
1013      !!
1014      !! ** Purpose : Read or write time-splitting arrays in restart file
1015      !!----------------------------------------------------------------------
1016      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1017      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1018      !
1019      !!----------------------------------------------------------------------
1020      !
1021      IF( TRIM(cdrw) == 'READ' ) THEN
1022         IF(nn_timing == 2)  CALL timing_start('iom_rstget')
1023         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
1024         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
1025         IF( .NOT.ln_bt_av ) THEN
1026            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
1027            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
1028            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
1029            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
1030            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
1031            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
1032         ENDIF
1033#if defined key_agrif
1034         ! Read time integrated fluxes
1035         IF ( .NOT.Agrif_Root() ) THEN
1036            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
1037            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
1038         ENDIF
1039#endif
1040      IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
1041      !
1042      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1043         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
1044         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
1045         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
1046         !
1047         IF (.NOT.ln_bt_av) THEN
1048            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
1049            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
1050            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
1051            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) )
1052            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) )
1053            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) )
1054         ENDIF
1055#if defined key_agrif
1056         ! Save time integrated fluxes
1057         IF ( .NOT.Agrif_Root() ) THEN
1058            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:) )
1059            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
1060         ENDIF
1061#endif
1062         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
1063      ENDIF
1064      !
1065   END SUBROUTINE ts_rst
1066
1067   SUBROUTINE dyn_spg_ts_init( kt )
1068      !!---------------------------------------------------------------------
1069      !!                   ***  ROUTINE dyn_spg_ts_init  ***
1070      !!
1071      !! ** Purpose : Set time splitting options
1072      !!----------------------------------------------------------------------
1073      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1074      !
1075      INTEGER  :: ji ,jj
1076      INTEGER  ::   ios                 ! Local integer output status for namelist read
1077      REAL(wp) :: zxr2, zyr2, zcmax
1078      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1079      !!
1080      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1081      &                  nn_baro, rn_bt_cmax, nn_bt_flt
1082      !!----------------------------------------------------------------------
1083      !
1084      REWIND( numnam_ref )              ! Namelist namsplit in reference namelist : time splitting parameters
1085      READ  ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1086901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1087
1088      REWIND( numnam_cfg )              ! Namelist namsplit in configuration namelist : time splitting parameters
1089      READ  ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1090902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1091      IF(lwm .AND. nprint > 2) WRITE ( numond, namsplit )
1092      !
1093      !         ! Max courant number for ext. grav. waves
1094      !
1095      CALL wrk_alloc( jpi, jpj, zcu )
1096      !
1097      IF (lk_vvl) THEN
1098         DO jj = 1, jpj
1099            DO ji =1, jpi
1100               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1101               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1102               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
1103            END DO
1104         END DO
1105      ELSE
1106         DO jj = 1, jpj
1107            DO ji =1, jpi
1108               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1109               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1110               zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )
1111            END DO
1112         END DO
1113      ENDIF
1114
1115      zcmax = MAXVAL(zcu(:,:))
1116      IF( lk_mpp )   CALL mpp_max( zcmax )
1117
1118      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1119      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1120     
1121      rdtbt = rdt / FLOAT(nn_baro)
1122      zcmax = zcmax * rdtbt
1123                     ! Print results
1124      IF(lwp) WRITE(numout,*)
1125      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1126      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1127      IF( ln_bt_nn_auto ) THEN
1128         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.true. Automatically set nn_baro '
1129         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
1130      ELSE
1131         IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1132      ENDIF
1133
1134      IF(ln_bt_av) THEN
1135         IF(lwp) WRITE(numout,*) '     ln_bt_av=.true.  => Time averaging over nn_baro time steps is on '
1136      ELSE
1137         IF(lwp) WRITE(numout,*) '     ln_bt_av=.false. => No time averaging of barotropic variables '
1138      ENDIF
1139      !
1140      !
1141      IF(ln_bt_fw) THEN
1142         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
1143      ELSE
1144         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
1145      ENDIF
1146      !
1147#if defined key_agrif
1148      ! Restrict the use of Agrif to the forward case only
1149      IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1150#endif
1151      !
1152      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
1153      SELECT CASE ( nn_bt_flt )
1154           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac'
1155           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
1156           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
1157           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1158      END SELECT
1159      !
1160      IF(lwp) WRITE(numout,*) ' '
1161      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
1162      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
1163      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1164      IF(lwp .AND. lflush) CALL flush(numout)
1165      !
1166      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1167         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1168      ENDIF
1169      IF ( zcmax>0.9_wp ) THEN
1170         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )         
1171      ENDIF
1172      !
1173      CALL wrk_dealloc( jpi, jpj, zcu )
1174      !
1175   END SUBROUTINE dyn_spg_ts_init
1176
1177#else
1178   !!---------------------------------------------------------------------------
1179   !!   Default case :   Empty module   No split explicit free surface
1180   !!---------------------------------------------------------------------------
1181CONTAINS
1182   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
1183      dyn_spg_ts_alloc = 0
1184   END FUNCTION dyn_spg_ts_alloc
1185   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
1186      INTEGER, INTENT(in) :: kt
1187      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
1188   END SUBROUTINE dyn_spg_ts
1189   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
1190      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1191      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
1192      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
1193   END SUBROUTINE ts_rst 
1194   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
1195      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
1196      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt
1197   END SUBROUTINE dyn_spg_ts_init
1198#endif
1199   
1200   !!======================================================================
1201END MODULE dynspg_ts
1202
1203
1204
Note: See TracBrowser for help on using the repository browser.