New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4438

Last change on this file since 4438 was 4433, checked in by trackstand2, 10 years ago

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