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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk

  • Property svn:keywords set to Id
File size: 42.3 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
4   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
5   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
6   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
7   !!              -   ! 2008-01  (R. Benshila)  change averaging method
8   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
10   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
11   !!---------------------------------------------------------------------
12#if defined key_dynspg_ts   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_dynspg_ts'         free surface cst volume with time splitting
15   !!----------------------------------------------------------------------
16   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
17   !!                 splitting scheme and add to the general trend
18   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE sbc_oce         ! surface boundary condition: ocean
23   USE dynspg_oce      ! surface pressure gradient variables
24   USE phycst          ! physical constants
25   USE domvvl          ! variable volume
26   USE zdfbfr          ! bottom friction
27   USE dynvor          ! vorticity term
28   USE obc_oce         ! Lateral open boundary condition
29   USE obc_par         ! open boundary condition parameters
30   USE obcdta          ! open boundary condition data     
31   USE obcfla          ! Flather open boundary condition 
32   USE bdy_par         ! for lk_bdy
33   USE bdy_oce         ! Lateral open boundary condition
34   USE bdydta          ! open boundary condition data     
35   USE bdydyn2d        ! open boundary conditions on barotropic variables
36   USE sbctide
37   USE updtide
38   USE lib_mpp         ! distributed memory computing library
39   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
40   USE prtctl          ! Print control
41   USE in_out_manager  ! I/O manager
42   USE iom             ! IOM library
43   USE restart         ! only for lrst_oce
44   USE zdf_oce         ! Vertical diffusion
45   USE wrk_nemo        ! Memory Allocation
46   USE timing          ! Timing
47
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC dyn_spg_ts        ! routine called by step.F90
53   PUBLIC ts_rst            ! routine called by istate.F90
54   PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90
55
56
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter
58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme)
59
60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b   ! now    averaged velocity
61   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b   ! before averaged velocity
62
63   !! * Substitutions
64#  include "domzgr_substitute.h90"
65#  include "vectopt_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
68   !! $Id$
69   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   INTEGER FUNCTION dyn_spg_ts_alloc()
74      !!----------------------------------------------------------------------
75      !!                  ***  routine dyn_spg_ts_alloc  ***
76      !!----------------------------------------------------------------------
77      ALLOCATE( ftnw  (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) ,     &
78         &      ftsw  (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc )
79         !
80      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc )
81      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
82      !
83   END FUNCTION dyn_spg_ts_alloc
84
85
86   SUBROUTINE dyn_spg_ts( kt )
87      !!----------------------------------------------------------------------
88      !!                  ***  routine dyn_spg_ts  ***
89      !!
90      !! ** Purpose :   Compute the now trend due to the surface pressure
91      !!      gradient in case of free surface formulation with time-splitting.
92      !!      Add it to the general trend of momentum equation.
93      !!
94      !! ** Method  :   Free surface formulation with time-splitting
95      !!      -1- Save the vertically integrated trend. This general trend is
96      !!          held constant over the barotropic integration.
97      !!          The Coriolis force is removed from the general trend as the
98      !!          surface gradient and the Coriolis force are updated within
99      !!          the barotropic integration.
100      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
101      !!          barotropic velocity (ua_e and va_e) through barotropic
102      !!          momentum and continuity integration. Barotropic former
103      !!          variables are time averaging over the full barotropic cycle
104      !!          (= 2 * baroclinic time step) and saved in uX_b
105      !!          and vX_b (X specifying after, now or before).
106      !!      -3- The new general trend becomes :
107      !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b )
108      !!
109      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
110      !!
111      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
112      !!---------------------------------------------------------------------
113      !
114      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
115      !
116      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
117      INTEGER  ::   icycle           ! local scalar
118      INTEGER  ::   ikbu, ikbv       ! local scalar
119      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars
120      REAL(wp) ::   z1_8, zx1, zy1                            !   -      -
121      REAL(wp) ::   z1_4, zx2, zy2                            !   -      -
122      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      -
123      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      -
124      REAL(wp) ::   ua_btm, va_btm                            !   -      -
125      !
126      REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv 
127      REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 
128      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum
129      !!----------------------------------------------------------------------
130      !
131      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts')
132      !
133      CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     )
134      CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
135      CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
136      !
137      IF( kt == nit000 ) THEN             !* initialisation
138         !
139         IF(lwp) WRITE(numout,*)
140         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
141         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
142         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
143         !
144         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b 
145         !
146         ua_e  (:,:) = un_b (:,:)
147         va_e  (:,:) = vn_b (:,:)
148         hu_e  (:,:) = hu   (:,:)
149         hv_e  (:,:) = hv   (:,:)
150         hur_e (:,:) = hur  (:,:)
151         hvr_e (:,:) = hvr  (:,:)
152         IF( ln_dynvor_een ) THEN
153            ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
154            DO jj = 2, jpj
155               DO ji = fs_2, jpi   ! vector opt.
156                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp
157                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp
158                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp
159                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp
160               END DO
161            END DO
162         ENDIF
163         !
164      ENDIF
165
166      !                                                     !* Local constant initialization
167      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step
168      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic
169                                                                        ! time step (euler timestep)
170      z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates
171      z1_4     = 0.25_wp       
172      zraur    = 1._wp / rau0                               ! 1 / volumic mass
173      !
174      zhdiv(:,:) = 0._wp                                    ! barotropic divergence
175      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F)
176      zv_sld = 0._wp   ;   zv_asp = 0._wp
177
178      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction
179        z2dt_bf = rdt
180      ELSE
181        z2dt_bf = 2.0_wp * rdt
182      ENDIF
183
184      ! -----------------------------------------------------------------------------
185      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
186      ! -----------------------------------------------------------------------------
187      !     
188      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
189      !                                   ! --------------------------
190      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp
191      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
192      !
193      DO jk = 1, jpkm1
194#if defined key_vectopt_loop
195         DO jj = 1, 1         !Vector opt. => forced unrolling
196            DO ji = 1, jpij
197#else
198         DO jj = 1, jpj
199            DO ji = 1, jpi
200#endif
201               !                                                                              ! now trend
202               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
203               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
204               !                                                                              ! now velocity
205               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
206               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
207               !
208#if defined key_vvl
209               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk) 
210               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk)
211#else
212               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
213               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
214#endif
215            END DO
216         END DO
217      END DO
218
219      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
220      DO jk = 1, jpkm1                    ! --------------------------
221         DO jj = 2, jpjm1
222            DO ji = fs_2, fs_jpim1   ! vector opt.
223               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
224               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
225            END DO
226         END DO
227      END DO
228
229      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
230      !                                   ! ---------------------------====
231      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport
232      zwy(:,:) = zvn(:,:) * e1v(:,:)
233      !
234      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1   ! vector opt.
237               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
238               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
239               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
240               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
241               ! energy conserving formulation for planetary vorticity term
242               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
243               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
244            END DO
245         END DO
246         !
247      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
248         DO jj = 2, jpjm1
249            DO ji = fs_2, fs_jpim1   ! vector opt.
250               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
251               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
252               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
253               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
254            END DO
255         END DO
256         !
257      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
258         DO jj = 2, jpjm1
259            DO ji = fs_2, fs_jpim1   ! vector opt.
260               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
261                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
262               zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
263                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
264            END DO
265         END DO
266         !
267      ENDIF
268
269      !                                   !* Right-Hand-Side of the barotropic momentum equation
270      !                                   ! ----------------------------------------------------
271      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
272         DO jj = 2, jpjm1 
273            DO ji = fs_2, fs_jpim1   ! vector opt.
274               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
275                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
276               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
277                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
278            END DO
279         END DO
280      ENDIF
281
282      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
283         DO ji = fs_2, fs_jpim1
284             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
285             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
286          END DO
287      END DO
288
289                   
290      !                                             ! Remove barotropic contribution of bottom friction
291      !                                             ! from the barotropic transport trend
292      zcoef = -1._wp * z1_2dt_b
293
294      IF(ln_bfrimp) THEN
295      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient
296      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction
297        DO jj = 2, jpjm1         
298           DO ji = fs_2, fs_jpim1
299              ikbu = mbku(ji,jj)
300              ikbv = mbkv(ji,jj)
301              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu)
302              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv)
303
304              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
305              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
306           END DO
307        END DO
308
309      ELSE
310
311# if defined key_vectopt_loop
312        DO jj = 1, 1
313           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
314# else
315        DO jj = 2, jpjm1
316           DO ji = 2, jpim1
317# endif
318            ! Apply stability criteria for bottom friction
319            !RBbug for vvl and external mode we may need to use varying fse3
320            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
321              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
322              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
323           END DO
324        END DO
325
326        IF( lk_vvl ) THEN
327           DO jj = 2, jpjm1
328              DO ji = fs_2, fs_jpim1   ! vector opt.
329                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
330                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) )
331                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
332                    &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) )
333              END DO
334           END DO
335        ELSE
336           DO jj = 2, jpjm1
337              DO ji = fs_2, fs_jpim1   ! vector opt.
338                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
339                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
340              END DO
341           END DO
342        ENDIF
343      END IF    ! end (ln_bfrimp)
344
345                   
346      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
347      !                                   ! --------------------------
348      zua(:,:) = zua(:,:) * hur(:,:)
349      zva(:,:) = zva(:,:) * hvr(:,:)
350      !
351      IF( lk_vvl ) THEN
352         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
353         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
354      ELSE
355         ub_b(:,:) = ub_b(:,:) * hur(:,:)
356         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
357      ENDIF
358
359      ! -----------------------------------------------------------------------
360      !  Phase 2 : Integration of the barotropic equations with time splitting
361      ! -----------------------------------------------------------------------
362      !
363      !                                             ! ==================== !
364      !                                             !    Initialisations   !
365      !                                             ! ==================== !
366      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
367     
368      !                                ! Start from NOW field
369      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
370      hv_e   (:,:) = hv   (:,:) 
371      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
372      hvr_e  (:,:) = hvr  (:,:)
373!RBbug     zsshb_e(:,:) = sshn (:,:) 
374      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
375      sshn_e (:,:) = sshn (:,:)
376     
377      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
378      zvn_e  (:,:) = vn_b (:,:)
379      zub_e  (:,:) = un_b (:,:)
380      zvb_e  (:,:) = vn_b (:,:)
381
382      zu_sum  (:,:) = un_b (:,:)           ! summation
383      zv_sum  (:,:) = vn_b (:,:)
384      zssh_sum(:,:) = sshn (:,:)
385
386#if defined key_obc
387      ! set ssh corrections to 0
388      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
389      IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp
390      IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp
391      IF( lp_obc_south )   sshfos_b(:,:) = 0._wp
392      IF( lp_obc_north )   sshfon_b(:,:) = 0._wp
393#endif
394
395      !                                             ! ==================== !
396      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
397         !                                          ! ==================== !
398         z2dt_e = 2. * ( rdt / nn_baro )
399         IF( jn == 1 )   z2dt_e = rdt / nn_baro
400
401         !                                                !* Update the forcing (BDY and tides)
402         !                                                !  ------------------
403         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
404         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
405         IF ( ln_tide_pot ) CALL upd_tide( kt, jn )
406
407         !                                                !* after ssh_e
408         !                                                !  -----------
409         DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports
410            DO ji = fs_2, fs_jpim1   ! vector opt.
411               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
412                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
413                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
414                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
415            END DO
416         END DO
417         !
418#if defined key_obc
419         !                                                     ! OBC : zhdiv must be zero behind the open boundary
420!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
421         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east
422         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west
423         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north
424         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south
425#endif
426#if defined key_bdy
427         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
428#endif
429         !
430         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
431            DO ji = fs_2, fs_jpim1   ! vector opt.
432               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
433            END DO
434         END DO
435
436         !                                                !* after barotropic velocities (vorticity scheme dependent)
437         !                                                !  --------------------------- 
438         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport
439         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
440         !
441         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
442            DO jj = 2, jpjm1
443               DO ji = fs_2, fs_jpim1   ! vector opt.
444                  ! surface pressure gradient
445                  IF( lk_vvl) THEN
446                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
447                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
448                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
449                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
450                  ELSE
451                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
452                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
453                  ENDIF
454                  ! add tidal astronomical forcing
455                  IF ( ln_tide_pot ) THEN
456                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
457                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
458                  ENDIF
459                  ! energy conserving formulation for planetary vorticity term
460                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
461                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
462                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
463                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
464                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
465                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
466                  ! after velocities with implicit bottom friction
467
468                  IF( ln_bfrimp ) THEN      ! implicit bottom friction
469                     !   A new method to implement the implicit bottom friction.
470                     !   H. Liu
471                     !   Sept 2011
472                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
473                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
474                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
475                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
476                     !
477                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
478                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
479                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
480                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
481                     !
482                  ELSE
483                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
484                      &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
485                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
486                      &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
487                  ENDIF
488               END DO
489            END DO
490            !
491         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
492            DO jj = 2, jpjm1
493               DO ji = fs_2, fs_jpim1   ! vector opt.
494                   ! surface pressure gradient
495                  IF( lk_vvl) THEN
496                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
497                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
498                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
499                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
500                  ELSE
501                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
502                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
503                  ENDIF
504                  ! add tidal astronomical forcing
505                  IF ( ln_tide_pot ) THEN
506                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
507                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
508                  ENDIF
509                  ! enstrophy conserving formulation for planetary vorticity term
510                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
511                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
512                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
513                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
514                  ! after velocities with implicit bottom friction
515                  IF( ln_bfrimp ) THEN
516                     !   A new method to implement the implicit bottom friction.
517                     !   H. Liu
518                     !   Sept 2011
519                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
520                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
521                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
522                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
523                     !
524                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
525                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
526                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
527                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
528                     !
529                  ELSE
530                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
531                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
532                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
533                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
534                  ENDIF
535               END DO
536            END DO
537            !
538         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
539            DO jj = 2, jpjm1
540               DO ji = fs_2, fs_jpim1   ! vector opt.
541                  ! surface pressure gradient
542                  IF( lk_vvl) THEN
543                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
544                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
545                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
546                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
547                  ELSE
548                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
549                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
550                  ENDIF
551                  ! add tidal astronomical forcing
552                  IF ( ln_tide_pot ) THEN
553                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
554                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
555                  ENDIF
556                  ! energy/enstrophy conserving formulation for planetary vorticity term
557                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
558                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
559                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
560                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
561                  ! after velocities with implicit bottom friction
562                  IF( ln_bfrimp ) THEN
563                     !   A new method to implement the implicit bottom friction.
564                     !   H. Liu
565                     !   Sept 2011
566                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            &
567                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            &
568                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
569                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
570                     !
571                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            &
572                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            &
573                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
574                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
575                     !
576                  ELSE
577                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   &
578                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
579                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   &
580                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
581                  ENDIF
582               END DO
583            END DO
584            !
585         ENDIF
586         !                                                     !* domain lateral boundary
587         !                                                     !  -----------------------
588
589                                                               ! OBC open boundaries
590         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
591
592                                                               ! BDY open boundaries
593#if defined key_bdy
594         pssh => sshn_e
595         phur => hur_e
596         phvr => hvr_e
597         pu2d => ua_e
598         pv2d => va_e
599
600         IF( lk_bdy )   CALL bdy_dyn2d( kt ) 
601#endif
602
603         !
604         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
605         CALL lbc_lnk( va_e  , 'V', -1. )
606         CALL lbc_lnk( ssha_e, 'T',  1. )
607
608         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
609         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
610         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
611
612         !                                                !* Time filter and swap
613         !                                                !  --------------------
614         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
615            zsshb_e(:,:) = sshn_e(:,:)
616            zub_e  (:,:) = zun_e (:,:)
617            zvb_e  (:,:) = zvn_e (:,:)
618            sshn_e (:,:) = ssha_e(:,:)
619            zun_e  (:,:) = ua_e  (:,:)
620            zvn_e  (:,:) = va_e  (:,:)
621         ELSE                                                   ! Swap + Filter
622            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
623            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
624            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
625            sshn_e (:,:) = ssha_e(:,:)
626            zun_e  (:,:) = ua_e  (:,:)
627            zvn_e  (:,:) = va_e  (:,:)
628         ENDIF
629
630         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
631            !                                             !  ------------------
632            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
633               DO ji = 1, fs_jpim1   ! Vector opt.
634                  zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
635                     &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
636                     &                     +   e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
637                  zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
638                     &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
639                     &                     +   e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
640               END DO
641            END DO
642            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
643            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
644            !
645            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
646            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
647            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
648            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
649            !
650         ENDIF
651         !                                                 ! ==================== !
652      END DO                                               !        end loop      !
653      !                                                    ! ==================== !
654
655#if defined key_obc
656      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
657      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
658      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
659      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
660#endif
661
662      ! -----------------------------------------------------------------------------
663      ! Phase 3. update the general trend with the barotropic trend
664      ! -----------------------------------------------------------------------------
665      !
666      !                                   !* Time average ==> after barotropic u, v, ssh
667      zcoef =  1._wp / ( 2 * nn_baro  + 1 ) 
668      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
669      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
670      !
671      !                                   !* update the general momentum trend
672      DO jk=1,jpkm1
673         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b
674         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b
675      END DO
676      un_b  (:,:) =  zu_sum(:,:) 
677      vn_b  (:,:) =  zv_sum(:,:) 
678      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
679      !
680      !                                   !* write time-spliting arrays in the restart
681      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
682      !
683      CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     )
684      CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   )
685      CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
686      !
687      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
688      !
689   END SUBROUTINE dyn_spg_ts
690
691
692   SUBROUTINE ts_rst( kt, cdrw )
693      !!---------------------------------------------------------------------
694      !!                   ***  ROUTINE ts_rst  ***
695      !!
696      !! ** Purpose : Read or write time-splitting arrays in restart file
697      !!----------------------------------------------------------------------
698      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
699      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
700      !
701      INTEGER ::  ji, jk        ! dummy loop indices
702      !!----------------------------------------------------------------------
703      !
704      IF( TRIM(cdrw) == 'READ' ) THEN
705         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
706            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
707            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
708         ELSE
709            un_b (:,:) = 0._wp
710            vn_b (:,:) = 0._wp
711            ! vertical sum
712            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
713               DO jk = 1, jpkm1
714                  DO ji = 1, jpij
715                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
716                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
717                  END DO
718               END DO
719            ELSE                             ! No  vector opt.
720               DO jk = 1, jpkm1
721                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
722                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
723               END DO
724            ENDIF
725            un_b (:,:) = un_b(:,:) * hur(:,:)
726            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
727         ENDIF
728
729         ! Vertically integrated velocity (before)
730         IF (neuler/=0) THEN
731            ub_b (:,:) = 0._wp
732            vb_b (:,:) = 0._wp
733
734            ! vertical sum
735            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
736               DO jk = 1, jpkm1
737                  DO ji = 1, jpij
738                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
739                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
740                  END DO
741               END DO
742            ELSE                             ! No  vector opt.
743               DO jk = 1, jpkm1
744                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
745                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
746               END DO
747            ENDIF
748
749            IF( lk_vvl ) THEN
750               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
751               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
752            ELSE
753               ub_b(:,:) = ub_b(:,:) * hur(:,:)
754               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
755            ENDIF
756         ELSE                                 ! neuler==0
757            ub_b (:,:) = un_b (:,:)
758            vb_b (:,:) = vn_b (:,:)
759         ENDIF
760
761         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
762            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh
763         ELSE
764            sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
765         ENDIF
766      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
767         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
768         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
769         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !
770      ENDIF
771      !
772   END SUBROUTINE ts_rst
773
774#else
775   !!----------------------------------------------------------------------
776   !!   Default case :   Empty module   No standart free surface cst volume
777   !!----------------------------------------------------------------------
778CONTAINS
779   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
780      dyn_spg_ts_alloc = 0
781   END FUNCTION dyn_spg_ts_alloc
782   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
783      INTEGER, INTENT(in) :: kt
784      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
785   END SUBROUTINE dyn_spg_ts
786   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
787      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
788      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
789      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
790   END SUBROUTINE ts_rst   
791#endif
792   
793   !!======================================================================
794END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.