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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 8 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

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