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

Last change on this file since 746 was 746, checked in by smasson, 16 years ago

implement ldstop in iom_varid, see ticket:21

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