source: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 788

Last change on this file since 788 was 788, checked in by rblod, 13 years ago

Change averaging method in dynspg_ts, see ticket #48

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