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

source: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 1405

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

VVL and time stplitting :

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