source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3876

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

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

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