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

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

source: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 1438

Last change on this file since 1438 was 1438, checked in by rblod, 15 years ago

Merge VVL branch with the trunk (act II), see ticket #429
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 30.4 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !! History :   1.0  !  04-12  (L. Bessieres, G. Madec)  Original code
4   !!              -   !  05-11  (V. Garnier, G. Madec)  optimization
5   !!              -   !  06-08  (S. Masson)  distributed restart using iom
6   !!             2.0  !  07-07  (D. Storkey) calls to BDY routines
7   !!              -   !  08-01  (R. Benshila)  change averaging method
8  !!---------------------------------------------------------------------
9#if defined key_dynspg_ts   ||   defined key_esopa
10   !!----------------------------------------------------------------------
11   !!   'key_dynspg_ts'         free surface cst volume with time splitting
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
15   !!                 splitting scheme and add to the general trend
16   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE sbc_oce         ! surface boundary condition: ocean
21   USE dynspg_oce      ! surface pressure gradient variables
22   USE phycst          ! physical constants
23   USE domvvl          ! variable volume
24   USE obcdta          ! open boundary condition data     
25   USE obcfla          ! Flather open boundary condition 
26   USE dynvor          ! vorticity term
27   USE obc_oce         ! Lateral open boundary condition
28   USE obc_par         ! open boundary condition parameters
29   USE bdy_oce         ! unstructured open boundaries
30   USE bdy_par         ! unstructured open boundaries
31   USE bdydta          ! unstructured open boundaries
32   USE bdydyn          ! unstructured open boundaries
33   USE bdytides        ! tidal forcing at unstructured open boundaries.
34   USE lib_mpp         ! distributed memory computing library
35   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
36   USE prtctl          ! Print control
37   USE in_out_manager  ! I/O manager
38   USE iom
39   USE restart         ! only for lrst_oce
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC dyn_spg_ts  ! routine called by step.F90
45   PUBLIC ts_rst      ! routine called by istate.F90
46
47   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter
48   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52#  include "vectopt_loop_substitute.h90"
53   !!-------------------------------------------------------------------------
54   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
55   !! $Id$
56   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
57   !!-------------------------------------------------------------------------
58
59CONTAINS
60
61   SUBROUTINE dyn_spg_ts( kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  routine dyn_spg_ts  ***
64      !!
65      !! ** Purpose :   Compute the now trend due to the surface pressure
66      !!      gradient in case of free surface formulation with time-splitting.
67      !!      Add it to the general trend of momentum equation.
68      !!      Compute the free surface.
69      !!
70      !! ** Method  :   Free surface formulation with time-splitting
71      !!      -1- Save the vertically integrated trend. This general trend is
72      !!          held constant over the barotropic integration.
73      !!          The Coriolis force is removed from the general trend as the
74      !!          surface gradient and the Coriolis force are updated within
75      !!          the barotropic integration.
76      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
77      !!          barotropic transports (ua_e and va_e) through barotropic
78      !!          momentum and continuity integration. Barotropic former
79      !!          variables are time averaging over the full barotropic cycle
80      !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b
81      !!          and zvX_b (X specifying after, now or before).
82      !!      -3- The new general trend becomes :
83      !!          ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H
84      !!
85      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
86      !!
87      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
88      !!---------------------------------------------------------------------
89      INTEGER, INTENT( in )  ::   kt           ! ocean time-step index
90      !!
91      INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices
92      INTEGER  ::  icycle                      ! temporary scalar
93      REAL(wp) ::                           &
94         zraur, zcoef, z2dt_e, z2dt_b, zfac25,   &  ! temporary scalars
95         zfact1, zspgu, zcubt, zx1, zy1,    &  !     "        "
96         zfact2, zspgv, zcvbt, zx2, zy2        !     "        "
97      REAL(wp), DIMENSION(jpi,jpj) ::       &
98         zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays
99         zua, zva, zub, zvb,                &  !     "        "
100         zua_e, zva_e, zssha_e,             &  !     "        "
101         zub_e, zvb_e, zsshb_e,             &  !     "        "
102         zu_sum, zv_sum, zssh_sum
103      !! Variable volume
104      REAL(wp), DIMENSION(jpi,jpj)     ::   &
105         zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace
106      !!----------------------------------------------------------------------
107
108      ! Arrays initialization
109      ! ---------------------
110      zhdiv(:,:) = 0.e0
111
112
113      IF( kt == nit000 ) THEN
114         !
115         IF(lwp) WRITE(numout,*)
116         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
117         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
118         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
119
120         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields:
121         !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
122
123         zssha_e(:,:) = sshn(:,:)
124         zua_e  (:,:) = un_e(:,:)
125         zva_e  (:,:) = vn_e(:,:)
126         hu_e   (:,:) = hu  (:,:)
127         hv_e   (:,:) = hv  (:,:)
128         IF( ln_dynvor_een ) THEN
129            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
130            DO jj = 2, jpj
131               DO ji = fs_2, jpi   ! vector opt.
132                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
133                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
134                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
135                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
136               END DO
137            END DO
138         ENDIF
139         !
140      ENDIF
141
142      ! Substract the surface pressure gradient from the momentum
143      ! ---------------------------------------------------------
144      IF( lk_vvl ) THEN    ! Variable volume
145
146         ! 0. Local initialization
147         ! -----------------------
148         zspgu_1(:,:) = 0.e0
149         zspgv_1(:,:) = 0.e0
150
151         ! 1. Surface pressure gradient (now)
152         ! ----------------------------
153         DO jj = 2, jpjm1
154            DO ji = fs_2, fs_jpim1   ! vector opt.
155               zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  ) &
156                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e1u(ji,jj)
157               zspgv_1(ji,jj) = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1) &
158                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e2v(ji,jj)
159            END DO
160         END DO
161
162         ! 2. Substract the surface pressure trend from the general trend
163         ! ------------------------------------------------------
164         DO jk = 1, jpkm1
165            DO jj = 2, jpjm1
166               DO ji = fs_2, fs_jpim1   ! vector opt.
167                  ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj)
168                  va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj)
169               END DO
170            END DO
171         END DO
172
173      ENDIF
174
175      ! Local constant initialization
176      ! --------------------------------
177      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
178      IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
179      zfact1 = 0.5 * 0.25                                   ! coefficient for vorticity estimates
180      zfact2 = 0.5 * 0.5
181      zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water
182
183      ! -----------------------------------------------------------------------------
184      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
185      ! -----------------------------------------------------------------------------
186
187      ! Vertically integrated quantities
188      ! --------------------------------
189      zua(:,:) = 0.e0
190      zva(:,:) = 0.e0
191      zub(:,:) = 0.e0
192      zvb(:,:) = 0.e0
193      zwx(:,:) = 0.e0
194      zwy(:,:) = 0.e0
195
196      ! vertical sum
197      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
198         DO jk = 1, jpkm1
199            DO ji = 1, jpij
200               !                                                           ! Vertically integrated momentum trends
201               zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk)
202               zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
203               !                                                           ! Vertically integrated transports (before)
204               zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
205               zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
206               !                                                           ! Planetary vorticity transport fluxes (now)
207               zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
208               zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk)
209            END DO
210         END DO
211      ELSE                             ! No  vector opt.
212         DO jk = 1, jpkm1
213            !                                                           ! Vertically integrated momentum trends
214            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
215            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
216            !                                                           ! Vertically integrated transports (before)
217            zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
218            zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
219            !                                                           ! Planetary vorticity (now)
220            zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
221            zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
222         END DO
223      ENDIF
224
225      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
226         DO jj = 2, jpjm1
227            DO ji = fs_2, fs_jpim1   ! vector opt.
228               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
229               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
230               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
231               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
232               ! energy conserving formulation for planetary vorticity term
233               zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
234               zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
235            END DO
236         END DO
237         !
238      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
239         DO jj = 2, jpjm1
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
242                              + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
243               zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
244                              + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
245               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
246               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
247            END DO
248         END DO
249         !
250      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
251         zfac25 = 0.25
252         DO jj = 2, jpjm1
253            DO ji = fs_2, fs_jpim1   ! vector opt.
254               zcu(ji,jj) = + zfac25 / e1u(ji,jj)   &
255                  &       * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
256                  &          + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
257               zcv(ji,jj) = - zfac25 / e2v(ji,jj)   &
258                  &       * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
259                  &          + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
260            END DO
261         END DO
262         !
263      ENDIF
264
265
266      ! Remove barotropic trend from general momentum trend
267      ! ---------------------------------------------------
268      DO jk = 1 , jpkm1
269         DO jj = 2, jpjm1
270            DO ji = fs_2, fs_jpim1   ! vector opt.
271               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
272               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
273            END DO
274         END DO
275      END DO
276
277      ! Remove coriolis term from barotropic trend
278      ! ------------------------------------------
279      DO jj = 2, jpjm1
280         DO ji = fs_2, fs_jpim1
281            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
282            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
283         END DO
284      END DO
285
286      ! -----------------------------------------------------------------------
287      !  Phase 2 : Integration of the barotropic equations with time splitting
288      ! -----------------------------------------------------------------------
289
290      ! Initialisations
291      !----------------
292      ! Number of iteration of the barotropic loop
293      icycle = 2  * nn_baro + 1
294
295      ! variables for the barotropic equations
296      zu_sum  (:,:) = 0.e0
297      zv_sum  (:,:) = 0.e0
298      zssh_sum(:,:) = 0.e0
299      hu_e    (:,:) = hu    (:,:)      ! (barotropic) ocean depth at u-point
300      hv_e    (:,:) = hv    (:,:)      ! (barotropic) ocean depth at v-point
301      zsshb_e (:,:) = sshn_e(:,:)      ! (barotropic) sea surface height (before and now)
302      ! vertical sum
303      un_e  (:,:) = 0.e0
304      vn_e  (:,:) = 0.e0
305      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
306         DO jk = 1, jpkm1
307            DO ji = 1, jpij
308               un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
309               vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
310            END DO
311         END DO
312      ELSE                             ! No  vector opt.
313         DO jk = 1, jpkm1
314            un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)
315            vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
316         END DO
317      ENDIF
318      zub_e (:,:) = un_e(:,:)
319      zvb_e (:,:) = vn_e(:,:)
320
321
322      ! set ssh corrections to 0
323      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
324#if defined key_obc
325      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
326      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
327      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
328      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
329#endif
330
331      ! Barotropic integration over 2 baroclinic time steps
332      ! ---------------------------------------------------
333
334      !                                                    ! ==================== !
335      DO jit = 1, icycle                                   !  sub-time-step loop  !
336         !                                                 ! ==================== !
337         z2dt_e = 2. * ( rdt / nn_baro )
338         IF ( jit == 1 )   z2dt_e = rdt / nn_baro
339
340         ! Time interpolation of open boundary condition data
341         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
342         IF( lk_bdy .OR. ln_bdy_tides)   CALL bdy_dta_bt( kt, jit+1 )
343
344         ! Horizontal divergence of barotropic transports
345         !--------------------------------------------------
346         zhdiv(:,:) = 0.e0
347         DO jj = 2, jpjm1
348            DO ji = fs_2, fs_jpim1   ! vector opt.
349               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * un_e(ji  ,jj)              &
350                  &            -e2u(ji-1,jj  ) * un_e(ji-1,jj)              &
351                  &            +e1v(ji  ,jj  ) * vn_e(ji  ,jj)              &
352                  &            -e1v(ji  ,jj-1) * vn_e(ji  ,jj-1) )          &
353                  &           / (e1t(ji,jj)*e2t(ji,jj))
354            END DO
355         END DO
356
357#if defined key_obc
358         ! open boundaries (div must be zero behind the open boundary)
359         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
360         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0      ! east
361         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0      ! west
362         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
363         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0      ! south
364#endif
365
366#if defined key_bdy
367         DO jj = 1, jpj
368            DO ji = 1, jpi
369               zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
370            END DO
371         END DO
372#endif
373
374         ! Sea surface height from the barotropic system
375         !----------------------------------------------
376         DO jj = 2, jpjm1
377            DO ji = fs_2, fs_jpim1   ! vector opt.
378               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
379                  &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
380            END DO
381         END DO
382
383         ! evolution of the barotropic transport ( following the vorticity scheme used)
384         ! ----------------------------------------------------------------------------
385         zwx(:,:) = e2u(:,:) * un_e(:,:)
386         zwy(:,:) = e1v(:,:) * vn_e(:,:)
387
388         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
389            DO jj = 2, jpjm1
390               DO ji = fs_2, fs_jpim1   ! vector opt.
391                  ! surface pressure gradient
392                  IF( lk_vvl) THEN
393                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) &
394                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj)
395                     zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) &
396                        &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj)
397                  ELSE
398                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
399                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
400                  ENDIF
401                  ! energy conserving formulation for planetary vorticity term
402                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
403                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
404                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
405                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
406                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
407                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
408                  ! after transports
409                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
410                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
411               END DO
412            END DO
413            !
414         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
415            DO jj = 2, jpjm1
416               DO ji = fs_2, fs_jpim1   ! vector opt.
417                  ! surface pressure gradient
418                  IF( lk_vvl) THEN
419                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) &
420                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj)
421                     zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) &
422                        &            - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj)
423                  ELSE
424                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
425                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
426                  ENDIF
427                  ! enstrophy conserving formulation for planetary vorticity term
428                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
429                     + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
430                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
431                     + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
432                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
433                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
434                  ! after transports
435                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
436                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
437               END DO
438            END DO
439            !
440         ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
441            zfac25 = 0.25
442            DO jj = 2, jpjm1
443               DO ji = fs_2, fs_jpim1   ! vector opt.
444                  ! surface pressure gradient
445                  IF( lk_vvl) THEN
446                     zspgu = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  ) &
447                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(ji,jj) / e1u(ji,jj)
448                     zspgv = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1) &
449                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hv(ji,jj) / e2v(ji,jj)
450                  ELSE
451                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
452                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
453                  ENDIF
454                  ! energy/enstrophy conserving formulation for planetary vorticity term
455                  zcubt = + zfac25 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
456                     &                             + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
457                  zcvbt = - zfac25 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
458                     &                             + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
459                  ! after transports
460                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
461                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
462               END DO
463            END DO
464            !
465         ENDIF
466
467         ! Flather's boundary condition for the barotropic loop :
468         !         - Update sea surface height on each open boundary
469         !         - Correct the barotropic transports
470         IF( lk_obc )   CALL obc_fla_ts
471         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla
472
473         ! ... Boundary conditions on ua_e, va_e, ssha_e
474         CALL lbc_lnk( zua_e  , 'U', -1. )
475         CALL lbc_lnk( zva_e  , 'V', -1. )
476         CALL lbc_lnk( zssha_e, 'T',  1. )
477
478         ! temporal sum
479         !-------------
480         zu_sum  (:,:) = zu_sum  (:,:) + zua_e  (:,:)
481         zv_sum  (:,:) = zv_sum  (:,:) + zva_e  (:,:) 
482         zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 
483
484         ! Time filter and swap of dynamics arrays
485         ! ---------------------------------------
486         IF( jit == 1 ) THEN   ! Euler (forward) time stepping
487            zsshb_e(:,:) = sshn_e (:,:)
488            zub_e  (:,:) = un_e   (:,:)
489            zvb_e  (:,:) = vn_e   (:,:)
490            sshn_e (:,:) = zssha_e(:,:)
491            un_e   (:,:) = zua_e  (:,:)
492            vn_e   (:,:) = zva_e  (:,:)
493         ELSE                                        ! Asselin filtering
494            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:)
495            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * un_e  (:,:)
496            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * vn_e  (:,:)
497            sshn_e (:,:) = zssha_e(:,:)
498            un_e   (:,:) = zua_e  (:,:)
499            vn_e   (:,:) = zva_e  (:,:)
500         ENDIF
501
502         IF( lk_vvl ) THEN ! Variable volume   ! needed for BDY ???
503
504            ! Sea surface elevation
505            ! ---------------------
506            DO jj = 1, jpjm1
507               DO ji = 1,jpim1
508
509                  ! Sea Surface Height at u-point before
510                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )   &
511                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
512                     &                + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
513
514                  ! Sea Surface Height at v-point before
515                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )   &
516                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
517                     &                + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
518
519               END DO
520            END DO
521
522            ! Boundaries conditions
523            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;  CALL lbc_lnk( zsshvn_e, 'V', 1. )
524
525            ! Ocean depth at U- and V-points
526            ! -------------------------------------------
527            hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:)
528            hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:)
529
530            !
531         ENDIF
532         !                                                 ! ==================== !
533      END DO                                               !        end loop      !
534      !                                                    ! ==================== !
535
536      ! Time average of after barotropic variables
537      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
538      un_e  (:,:) = zcoef *  zu_sum  (:,:) 
539      vn_e  (:,:) = zcoef *  zv_sum  (:,:) 
540      sshn_e(:,:) = zcoef *  zssh_sum(:,:) 
541
542#if defined key_obc
543      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
544      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
545      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
546      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
547#endif
548
549      ! -----------------------------------------------------------------------------
550      ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step)
551      ! -----------------------------------------------------------------------------
552
553
554
555      ! add time averaged barotropic coriolis and surface pressure gradient
556      ! terms to the general momentum trend
557      ! --------------------------------------------------------------------
558      DO jk=1,jpkm1
559         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b
560         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b
561      END DO
562
563      ! write filtered free surface arrays in restart file
564      ! --------------------------------------------------
565      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
566
567      !
568   END SUBROUTINE dyn_spg_ts
569
570
571   SUBROUTINE ts_rst( kt, cdrw )
572      !!---------------------------------------------------------------------
573      !!                   ***  ROUTINE ts_rst  ***
574      !!
575      !! ** Purpose : Read or write time-splitting arrays in restart file
576      !!----------------------------------------------------------------------
577      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
578      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
579      !
580      INTEGER ::  ji, jk        ! dummy loop indices
581      !!----------------------------------------------------------------------
582      !
583      IF( TRIM(cdrw) == 'READ' ) THEN
584         IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN
585            CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) )   ! free surface and
586            CALL iom_get( numror, jpdom_autoglo, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued
587            CALL iom_get( numror, jpdom_autoglo, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop
588         ELSE
589            sshn_e(:,:) = sshn(:,:)
590            un_e  (:,:) = 0.e0
591            vn_e  (:,:) = 0.e0
592            ! vertical sum
593            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
594               DO jk = 1, jpkm1
595                  DO ji = 1, jpij
596                     un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
597                     vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
598                  END DO
599               END DO
600            ELSE                             ! No  vector opt.
601               DO jk = 1, jpkm1
602                  un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)
603                  vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
604               END DO
605            ENDIF
606         ENDIF
607      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
608         CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) )   ! free surface and
609         CALL iom_rstput( kt, nitrst, numrow, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued
610         CALL iom_rstput( kt, nitrst, numrow, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop
611      ENDIF
612      !
613   END SUBROUTINE ts_rst
614
615#else
616   !!----------------------------------------------------------------------
617   !!   Default case :   Empty module   No standart free surface cst volume
618   !!----------------------------------------------------------------------
619CONTAINS
620   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
621      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
622   END SUBROUTINE dyn_spg_ts
623   SUBROUTINE ts_rst( kt, cdrw )     ! Empty routine
624      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
625      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
626      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
627   END SUBROUTINE ts_rst   
628#endif
629   
630   !!======================================================================
631END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.