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

Last change on this file since 553 was 553, checked in by opalod, 17 years ago

nemo_v1_bugfix_067: RB: mix-up between (zf)tne, (zf)tse, (zf)tnw, (zs)tsw for time-splitting
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 27.9 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   !!---------------------------------------------------------------------
10#if ( defined key_dynspg_ts && ! defined key_mpp_omp ) ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_dynspg_ts'         free surface cst volume with time splitting
13   !!   NOT 'key_mpp_omp'                          k-j-i loop (vector opt.)
14   !!----------------------------------------------------------------------
15   !!----------------------------------------------------------------------
16   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
17   !!                 splitting scheme and add to the general trend
18   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
19   !!----------------------------------------------------------------------
20   !! * Modules used
21   USE oce             ! ocean dynamics and tracers
22   USE dom_oce         ! ocean space and time domain
23   USE phycst          ! physical constants
24   USE ocesbc          ! ocean surface boundary condition
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 lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE prtctl          ! Print control
33   USE dynspg_oce      ! surface pressure gradient variables
34   USE in_out_manager  ! I/O manager
35   USE iom
36   USE restart         ! only for lrst_oce
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC dyn_spg_ts  ! routine called by step.F90
42
43   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne,   &  ! triad of coriolis parameter
44      &                             ftsw, ftse       ! (only used with een vorticity scheme)
45
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !!  OPA 9.0 , LOCEAN-IPSL (2005)
52   !! $Header$
53   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE dyn_spg_ts( kt )
59      !!----------------------------------------------------------------------
60      !!                  ***  routine dyn_spg_ts  ***
61      !!
62      !! ** Purpose :   Compute the now trend due to the surface pressure
63      !!      gradient in case of free surface formulation with time-splitting.
64      !!      Add it to the general trend of momentum equation.
65      !!      Compute the free surface.
66      !!
67      !! ** Method  :   Free surface formulation with time-splitting
68      !!      -1- Save the vertically integrated trend. This general trend is
69      !!          held constant over the barotropic integration.
70      !!          The Coriolis force is removed from the general trend as the
71      !!          surface gradient and the Coriolis force are updated within
72      !!          the barotropic integration.
73      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
74      !!          barotropic transports (ua_e and va_e) through barotropic
75      !!          momentum and continuity integration. Barotropic former
76      !!          variables are time averaging over the full barotropic cycle
77      !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b
78      !!          and zvX_b (X specifying after, now or before).
79      !!      -3- Update of sea surface height from time averaged barotropic
80      !!          variables.
81      !!        - apply lateral boundary conditions on sshn.
82      !!      -4- The new general trend becomes :
83      !!          ua = ua - sum_k(ua)/H + ( zua_b - 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      !! * Local declarations
92      INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices
93      INTEGER  ::  icycle                      ! temporary scalar
94      REAL(wp) ::                           &
95         zraur, zcoef, z2dt_e, z2dt_b, zfac25,   &  ! temporary scalars
96         zfact1, zspgu, zcubt, zx1, zy1,    &  !     "        "
97         zfact2, zspgv, zcvbt, zx2, zy2        !     "        "
98      REAL(wp), DIMENSION(jpi,jpj) ::       &
99         zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays
100         zua, zva, zub, zvb,                &  !     "        "
101         zssha_b, zua_b, zva_b,             &  !     "        "
102         zsshb_e, zub_e, zvb_e,             &  !     "        "
103         zun_e, zvn_e                          !     "        "
104      !!----------------------------------------------------------------------
105
106      ! Arrays initialization
107      ! ---------------------
108      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0
109      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0
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 = ', FLOOR( 2*rdt/rdtbt )
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         ssha_e(:,:) = sshn(:,:)
124         ua_e(:,:)   = un_b(:,:)
125         va_e(:,:)   = vn_b(:,:)
126
127         IF( ln_dynvor_een ) THEN
128            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
129            DO jj = 2, jpj
130               DO ji = fs_2, jpi   ! vector opt.
131                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
132                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
133                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
134                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
135               END DO
136            END DO
137         ENDIF
138         !
139      ENDIF
140
141      ! Local constant initialization
142      ! --------------------------------
143      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
144      IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
145      zfact1 = 0.5 * 0.25                                   ! coefficient for vorticity estimates
146      zfact2 = 0.5 * 0.5
147      zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water
148     
149      ! -----------------------------------------------------------------------------
150      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
151      ! -----------------------------------------------------------------------------
152
153      ! Vertically integrated quantities
154      ! --------------------------------
155      zua(:,:) = 0.e0
156      zva(:,:) = 0.e0
157      zub(:,:) = 0.e0
158      zvb(:,:) = 0.e0
159      zwx(:,:) = 0.e0
160      zwy(:,:) = 0.e0
161
162      ! vertical sum
163      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
164         DO jk = 1, jpkm1
165            DO ji = 1, jpij
166               !                                                           ! Vertically integrated momentum trends
167               zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk)
168               zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
169               !                                                           ! Vertically integrated transports (before)
170               zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk)
171               zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk)
172               !                                                           ! Planetary vorticity transport fluxes (now)
173               zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
174               zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk)
175            END DO
176         END DO
177      ELSE                             ! No  vector opt.
178         DO jk = 1, jpkm1
179            !                                                           ! Vertically integrated momentum trends
180            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
181            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
182            !                                                           ! Vertically integrated transports (before)
183            zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk)
184            zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk)
185            !                                                           ! Planetary vorticity (now)
186            zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
187            zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
188         END DO
189      ENDIF
190
191      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
192         DO jj = 2, jpjm1
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
195               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
196               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
197               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
198               ! energy conserving formulation for planetary vorticity term
199               zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
200               zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
201            END DO
202         END DO
203         !
204      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
205         DO jj = 2, jpjm1
206            DO ji = fs_2, fs_jpim1   ! vector opt.
207               zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
208                              + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
209               zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
210                              + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
211               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
212               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
213            END DO
214         END DO
215         !
216      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
217         zfac25 = 0.25
218         DO jj = 2, jpjm1
219            DO ji = fs_2, fs_jpim1   ! vector opt.
220               zcu(ji,jj) = + zfac25 / e1u(ji,jj)   &
221                  &       * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
222                  &          + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
223               zcv(ji,jj) = - zfac25 / e2v(ji,jj)   &
224                  &       * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
225                  &          + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
226            END DO
227         END DO
228         !
229      ENDIF
230
231
232      ! Remove barotropic trend from general momentum trend
233      ! ---------------------------------------------------
234      DO jk = 1 , jpkm1
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1   ! vector opt.
237               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
238               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
239            END DO
240         END DO
241      END DO
242
243      ! Remove coriolis term from barotropic trend
244      ! ------------------------------------------
245      DO jj = 2, jpjm1
246         DO ji = fs_2, fs_jpim1
247            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
248            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
249         END DO
250      END DO
251
252      ! -----------------------------------------------------------------------
253      !  Phase 2 : Integration of the barotropic equations with time splitting
254      ! -----------------------------------------------------------------------
255
256      ! Initialisations
257      !----------------
258      ! Number of iteration of the barotropic loop
259      icycle = FLOOR( z2dt_b / rdtbt )
260
261      ! variables for the barotropic equations
262      zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
263      sshn_e (:,:) = sshn_b(:,:)
264      zub_e  (:,:) = un_b  (:,:)       ! barotropic transports issued from the barotropic equations (before and now)
265      zvb_e  (:,:) = vn_b  (:,:)
266      zun_e  (:,:) = un_b  (:,:)
267      zvn_e  (:,:) = vn_b  (:,:)
268      zssha_b(:,:) = sshn  (:,:)       ! time averaged variables over all sub-timesteps
269      zua_b  (:,:) = un_b  (:,:)   
270      zva_b  (:,:) = vn_b  (:,:)
271
272      ! set ssh corrections to 0
273      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
274#if defined key_obc
275      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
276      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
277      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
278      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
279#endif
280
281      ! Barotropic integration over 2 baroclinic time steps
282      ! ---------------------------------------------------
283
284      !                                                    ! ==================== !
285      DO jit = 1, icycle                                   !  sub-time-step loop  !
286         !                                                 ! ==================== !
287         z2dt_e = 2. * rdtbt
288         IF ( jit == 1 )   z2dt_e = rdtbt
289
290         ! Time interpolation of open boundary condition data
291         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
292
293         ! Horizontal divergence of barotropic transports
294         !--------------------------------------------------
295         DO jj = 2, jpjm1
296            DO ji = fs_2, fs_jpim1   ! vector opt.
297               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              &
298                  &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              &
299                  &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              &
300                  &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          &
301                  &           / (e1t(ji,jj)*e2t(ji,jj))
302            END DO
303         END DO
304
305#if defined key_obc
306         ! open boundaries (div must be zero behind the open boundary)
307         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
308         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0      ! east
309         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0      ! west
310         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
311         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0      ! south
312#endif
313
314         ! Sea surface height from the barotropic system
315         !----------------------------------------------
316         DO jj = 2, jpjm1
317            DO ji = fs_2, fs_jpim1   ! vector opt.
318               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
319            &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
320            END DO
321         END DO
322
323         ! evolution of the barotropic transport ( following the vorticity scheme used)
324         ! ----------------------------------------------------------------------------
325         zwx(:,:) = e2u(:,:) * zun_e(:,:)
326         zwy(:,:) = e1v(:,:) * zvn_e(:,:)
327
328         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
329            DO jj = 2, jpjm1
330               DO ji = fs_2, fs_jpim1   ! vector opt.
331                  ! surface pressure gradient
332                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
333                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
334                  ! energy conserving formulation for planetary vorticity term
335                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
336                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
337                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
338                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
339                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
340                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
341                  ! after transports
342                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
343                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
344               END DO
345            END DO
346            !
347         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
348            DO jj = 2, jpjm1
349               DO ji = fs_2, fs_jpim1   ! vector opt.
350                  ! surface pressure gradient
351                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
352                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
353                  ! enstrophy conserving formulation for planetary vorticity term
354                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
355                                 + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
356                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
357                                 + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
358                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
359                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
360                  ! after transports
361                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
362                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
363               END DO
364            END DO
365            !
366         ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
367            zfac25 = 0.25
368            DO jj = 2, jpjm1
369               DO ji = fs_2, fs_jpim1   ! vector opt.
370                  ! surface pressure gradient
371                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
372                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
373                  ! energy/enstrophy conserving formulation for planetary vorticity term
374                  zcubt = + zfac25 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
375                     &                             + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
376                  zcvbt = - zfac25 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
377                     &                             + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
378                  ! after transports
379                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
380                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
381               END DO
382            END DO
383            !
384         ENDIF
385
386         ! Flather's boundary condition for the barotropic loop :
387         !         - Update sea surface height on each open boundary
388         !         - Correct the barotropic transports
389         IF( lk_obc )   CALL obc_fla_ts
390
391
392         ! ... Boundary conditions on ua_e, va_e, ssha_e
393         CALL lbc_lnk( ua_e  , 'U', -1. )
394         CALL lbc_lnk( va_e  , 'V', -1. )
395         CALL lbc_lnk( ssha_e, 'T',  1. )
396
397         ! temporal sum
398         !-------------
399         zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:)
400         zua_b  (:,:) = zua_b  (:,:) + ua_e  (:,:)
401         zva_b  (:,:) = zva_b  (:,:) + va_e  (:,:) 
402
403         ! Time filter and swap of dynamics arrays
404         ! ---------------------------------------
405         IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
406            zsshb_e(:,:) = sshn_e(:,:)
407            zub_e  (:,:) = zun_e (:,:)
408            zvb_e  (:,:) = zvn_e (:,:)
409            sshn_e (:,:) = ssha_e(:,:)
410            zun_e  (:,:) = ua_e  (:,:)
411            zvn_e  (:,:) = va_e  (:,:)
412         ELSE                                        ! Asselin filtering
413            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
414            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:)
415            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:)
416            sshn_e (:,:) = ssha_e(:,:)
417            zun_e  (:,:) = ua_e  (:,:)
418            zvn_e  (:,:) = va_e  (:,:)
419         ENDIF
420
421         !                                                 ! ==================== !
422      END DO                                               !        end loop      !
423      !                                                    ! ==================== !
424
425
426      ! Time average of after barotropic variables
427      zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
428      zssha_b(:,:) = zcoef * zssha_b(:,:) 
429      zua_b  (:,:) = zcoef *  zua_b (:,:) 
430      zva_b  (:,:) = zcoef *  zva_b (:,:) 
431#if defined key_obc
432         IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
433         IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
434         IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
435         IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
436#endif
437     
438
439      ! ---------------------------------------------------------------------------
440      ! Phase 3 : Update sea surface height from time averaged barotropic variables
441      ! ---------------------------------------------------------------------------
442
443 
444      ! Horizontal divergence of time averaged barotropic transports
445      !-------------------------------------------------------------
446      DO jj = 2, jpjm1
447         DO ji = fs_2, fs_jpim1   ! vector opt.
448            zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
449           &                +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
450           &             / ( e1t(ji,jj) * e2t(ji,jj) )
451         END DO
452      END DO
453
454#if defined key_obc
455      ! open boundaries (div must be zero behind the open boundary)
456      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
457      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
458      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
459      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
460      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
461#endif
462
463      ! sea surface height
464      !-------------------
465      sshb(:,:) = sshn(:,:)
466      sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
467
468      ! ... Boundary conditions on sshn
469      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
470
471
472      ! -----------------------------------------------------------------------------
473      ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
474      ! -----------------------------------------------------------------------------
475
476      ! Swap on time averaged barotropic variables
477      ! ------------------------------------------
478      sshb_b(:,:) = sshn_b (:,:)
479      sshn_b(:,:) = zssha_b(:,:)
480      un_b  (:,:) = zua_b  (:,:) 
481      vn_b  (:,:) = zva_b  (:,:) 
482   
483      ! add time averaged barotropic coriolis and surface pressure gradient
484      ! terms to the general momentum trend
485      ! --------------------------------------------------------------------
486      DO jk=1,jpkm1
487         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b
488         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b
489      END DO
490
491      ! write filtered free surface arrays in restart file
492      ! --------------------------------------------------
493      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
494
495      ! print sum trends (used for debugging)
496      IF( ln_ctl )     CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
497      !
498   END SUBROUTINE dyn_spg_ts
499
500
501   SUBROUTINE ts_rst( kt, cdrw )
502      !!---------------------------------------------------------------------
503      !!                   ***  ROUTINE ts_rst  ***
504      !!
505      !! ** Purpose : Read or write time-splitting arrays in restart file
506      !!----------------------------------------------------------------------
507      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
508      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
509      !
510      INTEGER ::  ji, jk        ! dummy loop indices
511      !!----------------------------------------------------------------------
512      !
513      IF( TRIM(cdrw) == 'READ' ) THEN
514         IF( iom_varid( numror, 'sshn' ) > 0 ) THEN
515            CALL iom_get( numror, jpdom_local, 'sshb'  , sshb(:,:)   )
516            CALL iom_get( numror, jpdom_local, 'sshn'  , sshn(:,:)   )
517            IF( neuler == 0 ) sshb(:,:) = sshn(:,:)
518         ELSE
519            sshb(:,:) = 0.e0
520            sshn(:,:) = 0.e0
521         ENDIF
522         IF( iom_varid( numror, 'sshn_b' ) > 0 ) THEN
523            CALL iom_get( numror, jpdom_local, 'sshb_b', sshb_b(:,:) )   ! free surface issued
524            CALL iom_get( numror, jpdom_local, 'sshn_b', sshn_b(:,:) )   ! from time-splitting loop
525            CALL iom_get( numror, jpdom_local, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
526            CALL iom_get( numror, jpdom_local, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
527            IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:)
528         ELSE
529            sshb_b(:,:) = sshb(:,:)
530            sshn_b(:,:) = sshn(:,:)
531            un_b  (:,:) = 0.e0
532            vn_b  (:,:) = 0.e0
533            ! vertical sum
534            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
535               DO jk = 1, jpkm1
536                  DO ji = 1, jpij
537                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
538                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
539                  END DO
540               END DO
541            ELSE                             ! No  vector opt.
542               DO jk = 1, jpkm1
543                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
544                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
545               END DO
546            ENDIF
547         ENDIF
548      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
549         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) )
550         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) )
551         CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) )   ! free surface issued
552         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) )   ! from barotropic loop
553         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
554         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
555      ENDIF
556      !
557   END SUBROUTINE ts_rst
558
559#else
560   !!----------------------------------------------------------------------
561   !!   Default case :   Empty module   No standart free surface cst volume
562   !!----------------------------------------------------------------------
563CONTAINS
564   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
565      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
566   END SUBROUTINE dyn_spg_ts
567#endif
568   
569   !!======================================================================
570END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.