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

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

nemo_v1_bugfix_079: CT : - add the creation of restart file (saving sshn & sshb) when using the explicit surface pressure option

  • light modifications in the way to read/write restart files in dynspg_* and associated _jki_ modules (key_mpp_omp removed)
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 28.0 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_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_dynspg_ts'         free surface cst volume with time splitting
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
16   !!                 splitting scheme and add to the general trend
17   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
18   !!----------------------------------------------------------------------
19   !! * Modules used
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE phycst          ! physical constants
23   USE ocesbc          ! ocean surface boundary condition
24   USE obcdta          ! open boundary condition data     
25   USE obcfla          ! Flather open boundary condition 
26   USE dynvor          ! vorticity term
27   USE obc_oce         ! Lateral open boundary condition
28   USE obc_par         ! open boundary condition parameters
29   USE lib_mpp         ! distributed memory computing library
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE prtctl          ! Print control
32   USE dynspg_oce      ! surface pressure gradient variables
33   USE in_out_manager  ! I/O manager
34   USE iom
35   USE restart         ! only for lrst_oce
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC dyn_spg_ts  ! routine called by step.F90
41   PUBLIC ts_rst      ! routine called by j-k-i subroutine
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            IF( nn_rstssh == 1 ) THEN 
520               sshb(:,:) = 0.e0
521               sshn(:,:) = 0.e0
522            ENDIF
523         ENDIF
524         IF( iom_varid( numror, 'sshn_b' ) > 0 ) THEN
525            CALL iom_get( numror, jpdom_local, 'sshb_b', sshb_b(:,:) )   ! free surface issued
526            CALL iom_get( numror, jpdom_local, 'sshn_b', sshn_b(:,:) )   ! from time-splitting loop
527            CALL iom_get( numror, jpdom_local, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
528            CALL iom_get( numror, jpdom_local, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
529            IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:)
530         ELSE
531            sshb_b(:,:) = sshb(:,:)
532            sshn_b(:,:) = sshn(:,:)
533            un_b  (:,:) = 0.e0
534            vn_b  (:,:) = 0.e0
535            ! vertical sum
536            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
537               DO jk = 1, jpkm1
538                  DO ji = 1, jpij
539                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
540                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
541                  END DO
542               END DO
543            ELSE                             ! No  vector opt.
544               DO jk = 1, jpkm1
545                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
546                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
547               END DO
548            ENDIF
549         ENDIF
550      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
551         CALL iom_rstput( kt, nitrst, numrow, 'sshb'  , sshb  (:,:) )
552         CALL iom_rstput( kt, nitrst, numrow, 'sshn'  , sshn  (:,:) )
553         CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) )   ! free surface issued
554         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) )   ! from barotropic loop
555         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
556         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
557      ENDIF
558      !
559   END SUBROUTINE ts_rst
560
561#else
562   !!----------------------------------------------------------------------
563   !!   Default case :   Empty module   No standart free surface cst volume
564   !!----------------------------------------------------------------------
565CONTAINS
566   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
567      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
568   END SUBROUTINE dyn_spg_ts
569#endif
570   
571   !!======================================================================
572END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.