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

source: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7171

Last change on this file since 7171 was 7171, checked in by emanuelaclementi, 7 years ago

#1643 Updates and bug fix on development branch 2015/dev_r5936_INGV1_WAVE

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