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/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 11277

Last change on this file since 11277 was 11277, checked in by kingr, 3 years ago

Merged Juan's changes for running AMM15 woth wave coupling.
Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results.
Few line spacing changes to allow merging with OBS br and trunk rvn 5518.

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