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

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

source: branches/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5427

Last change on this file since 5427 was 5427, checked in by deazer, 9 years ago

Added back in basic shelf seas diagnostics after removal of svn keywords.
Builds, extracts and mereges and runs as expected from working copy

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