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 @ 1241

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

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