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

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

Correction of major bug on bottom friction, see ticket #233

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 32.3 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
4   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
5   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
6   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
7   !!              -   ! 2008-01  (R. Benshila)  change averaging method
8   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9  !!---------------------------------------------------------------------
10#if defined key_dynspg_ts   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_dynspg_ts'         free surface cst volume with time splitting
13   !!----------------------------------------------------------------------
14   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
15   !!                 splitting scheme and add to the general trend
16   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE sbc_oce         ! surface boundary condition: ocean
21   USE dynspg_oce      ! surface pressure gradient variables
22   USE phycst          ! physical constants
23   USE domvvl          ! variable volume
24   USE zdfbfr          ! bottom friction
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 bdy_oce         ! unstructured open boundaries
31   USE bdy_par         ! unstructured open boundaries
32   USE bdydta          ! unstructured open boundaries
33   USE bdydyn          ! unstructured open boundaries
34   USE bdytides        ! tidal forcing at unstructured open boundaries.
35   USE lib_mpp         ! distributed memory computing library
36   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
37   USE prtctl          ! Print control
38   USE in_out_manager  ! I/O manager
39   USE iom
40   USE restart         ! only for lrst_oce
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC dyn_spg_ts  ! routine called by step.F90
46   PUBLIC ts_rst      ! routine called by istate.F90
47
48
49   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter
50   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
51
52   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! averaged velocity
53
54   !! * Substitutions
55#  include "domzgr_substitute.h90"
56#  include "vectopt_loop_substitute.h90"
57   !!-------------------------------------------------------------------------
58   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
59   !! $Id$
60   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
61   !!-------------------------------------------------------------------------
62
63CONTAINS
64
65   SUBROUTINE dyn_spg_ts( kt )
66      !!----------------------------------------------------------------------
67      !!                  ***  routine dyn_spg_ts  ***
68      !!
69      !! ** Purpose :   Compute the now trend due to the surface pressure
70      !!      gradient in case of free surface formulation with time-splitting.
71      !!      Add it to the general trend of momentum equation.
72      !!
73      !! ** Method  :   Free surface formulation with time-splitting
74      !!      -1- Save the vertically integrated trend. This general trend is
75      !!          held constant over the barotropic integration.
76      !!          The Coriolis force is removed from the general trend as the
77      !!          surface gradient and the Coriolis force are updated within
78      !!          the barotropic integration.
79      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
80      !!          barotropic velocity (ua_e and va_e) through barotropic
81      !!          momentum and continuity integration. Barotropic former
82      !!          variables are time averaging over the full barotropic cycle
83      !!          (= 2 * baroclinic time step) and saved in zuX_b
84      !!          and zvX_b (X specifying after, now or before).
85      !!      -3- The new general trend becomes :
86      !!          ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) )
87      !!
88      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
89      !!
90      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
91      !!---------------------------------------------------------------------
92      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
93      !!
94      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
95      INTEGER  ::   icycle           ! temporary scalar
96
97      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! temporary scalars
98      REAL(wp) ::   z1_8, zx1, zy1                   !    -         -
99      REAL(wp) ::   z1_4, zx2, zy2                   !     -         -
100      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !     -         -
101      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !     -         -
102      !!
103      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e
104      !!
105      REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace
106      !!
107      REAL(wp), DIMENSION(jpi,jpj) ::   zcu, zwx, zua, zun, zub   ! 2D workspace
108      REAL(wp), DIMENSION(jpi,jpj) ::   zcv, zwy, zva, zvn, zvb   !  -      -
109      REAL(wp), DIMENSION(jpi,jpj) ::   zun_e, zub_e, zu_sum      ! 2D workspace
110      REAL(wp), DIMENSION(jpi,jpj) ::   zvn_e, zvb_e, zv_sum      !  -      -
111      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_sum                  !  -      -
112      !!----------------------------------------------------------------------
113
114      IF( kt == nit000 ) THEN             !* initialisation
115         !
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
118         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
119         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
120         !
121         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields:
122         !                               ! un_b, vn_b 
123         !
124         ua_e  (:,:) = un_b (:,:)
125         va_e  (:,:) = vn_b (:,:)
126         hu_e  (:,:) = hu   (:,:)
127         hv_e  (:,:) = hv   (:,:)
128         hur_e (:,:) = hur  (:,:)
129         hvr_e (:,:) = hvr  (:,:)
130         IF( ln_dynvor_een ) THEN
131            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
132            DO jj = 2, jpj
133               DO ji = fs_2, jpi   ! vector opt.
134                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
135                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
136                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
137                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
138               END DO
139            END DO
140         ENDIF
141         !
142      ENDIF
143
144      !                                   !* Local constant initialization
145      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
146      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
147      z1_4 = 0.5 * 0.5
148      zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water
149      !
150      zhdiv(:,:) = 0.e0                                     ! barotropic divergence
151      zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F)
152      zv_sld = 0.e0   ;   zv_asp = 0.e0
153
154      ! -----------------------------------------------------------------------------
155      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
156      ! -----------------------------------------------------------------------------
157      !     
158      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
159      !                                   ! --------------------------
160      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   zub(:,:) = 0.e0
161      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   zvb(:,:) = 0.e0
162      !
163      DO jk = 1, jpkm1
164#if defined key_vectopt_loop
165         DO jj = 1, 1         !Vector opt. => forced unrolling
166            DO ji = 1, jpij
167#else
168         DO jj = 1, jpj
169            DO ji = 1, jpi
170#endif
171               !                                                                              ! now trend
172               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
173               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
174               !                                                                              ! now velocity
175               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
176               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
177               !                                                                              ! before velocity
178               zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 
179               zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
180            END DO
181         END DO
182      END DO
183
184      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
185      DO jk = 1, jpkm1                    ! --------------------------
186         DO jj = 2, jpjm1
187            DO ji = fs_2, fs_jpim1   ! vector opt.
188               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
189               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
190            END DO
191         END DO
192      END DO
193
194      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
195      !                                   ! ---------------------------====
196      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport
197      zwy(:,:) = zvn(:,:) * e1v(:,:)
198      !
199      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
200         DO jj = 2, jpjm1
201            DO ji = fs_2, fs_jpim1   ! vector opt.
202               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
203               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
204               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
205               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
206               ! energy conserving formulation for planetary vorticity term
207               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
208               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
209            END DO
210         END DO
211         !
212      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
213         DO jj = 2, jpjm1
214            DO ji = fs_2, fs_jpim1   ! vector opt.
215               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
216               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
217               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
218               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
219            END DO
220         END DO
221         !
222      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
223         DO jj = 2, jpjm1
224            DO ji = fs_2, fs_jpim1   ! vector opt.
225               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
226                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
227               zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
228                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
229            END DO
230         END DO
231         !
232      ENDIF
233
234      !                                   !* Right-Hand-Side of the barotropic momentum equation
235      !                                   ! ----------------------------------------------------
236      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
237         DO jj = 2, jpjm1 
238            DO ji = fs_2, fs_jpim1   ! vector opt.
239               zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    &
240                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
241               zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    &
242                  &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
243            END DO
244         END DO
245      ENDIF
246
247      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
248         DO ji = fs_2, fs_jpim1
249            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
250            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
251         END DO
252      END DO
253
254      !                                            ! Remove barotropic contribution of bottom friction
255                                                   ! from the barotropic transport trend
256      IF( lk_vvl ) THEN
257       DO jj = 2, jpjm1
258          DO ji = fs_2, fs_jpim1   ! vector opt.
259             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
260             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
261          END DO
262       END DO
263      ELSE
264       DO jj = 2, jpjm1
265          DO ji = fs_2, fs_jpim1   ! vector opt.
266             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) * hur(ji,jj)
267             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) * hvr(ji,jj)
268          END DO
269       END DO
270      ENDIF
271
272      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
273      !                                   ! --------------------------
274      zua(:,:) = zua(:,:) * hur(:,:)
275      zva(:,:) = zva(:,:) * hvr(:,:)
276      !
277      IF( lk_vvl ) THEN
278         zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
279         zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
280      ELSE
281         zub(:,:) = zub(:,:) * hur(:,:)
282         zvb(:,:) = zvb(:,:) * hvr(:,:)
283      ENDIF
284
285      ! -----------------------------------------------------------------------
286      !  Phase 2 : Integration of the barotropic equations with time splitting
287      ! -----------------------------------------------------------------------
288      !
289      !                                             ! ==================== !
290      !                                             !    Initialisations   !
291      !                                             ! ==================== !
292      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
293     
294      !                                ! Start from NOW field
295      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
296      hv_e   (:,:) = hv   (:,:) 
297      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
298      hvr_e  (:,:) = hvr  (:,:)
299!RBbug     zsshb_e(:,:) = sshn (:,:) 
300      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
301      sshn_e (:,:) = sshn (:,:)
302     
303      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
304      zvn_e  (:,:) = vn_b (:,:)
305      zub_e  (:,:) = un_b (:,:)
306      zvb_e  (:,:) = vn_b (:,:)
307
308      zu_sum  (:,:) = un_b (:,:)           ! summation
309      zv_sum  (:,:) = vn_b (:,:)
310      zssh_sum(:,:) = sshn (:,:)
311
312#if defined key_obc
313      ! set ssh corrections to 0
314      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
315      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
316      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
317      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
318      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
319#endif
320
321      !                                             ! ==================== !
322      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
323         !                                          ! ==================== !
324         z2dt_e = 2. * ( rdt / nn_baro )
325         IF( jn == 1 )   z2dt_e = rdt / nn_baro
326
327         !                                                !* Update the forcing (OBC, BDY and tides)
328         !                                                !  ------------------
329         IF( lk_obc                     )   CALL obc_dta_bt( kt, jn   )
330         IF( lk_bdy  .OR.  ln_bdy_tides )   CALL bdy_dta_bt( kt, jn+1 )
331
332         !                                                !* after ssh_e
333         !                                                !  -----------
334         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
335            DO ji = fs_2, fs_jpim1   ! vector opt.
336               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
337                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
338                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
339                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
340            END DO
341         END DO
342         !
343#if defined key_obc
344         !                                                     ! OBC : zhdiv must be zero behind the open boundary
345!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
346         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
347         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
348         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
349         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
350#endif
351#if defined key_bdy
352         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
353#endif
354         !
355         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
356            DO ji = fs_2, fs_jpim1   ! vector opt.
357               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
358            END DO
359         END DO
360
361         !                                                !* after barotropic velocities (vorticity scheme dependent)
362         !                                                !  --------------------------- 
363         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
364         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
365         !
366         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
367            DO jj = 2, jpjm1
368               DO ji = fs_2, fs_jpim1   ! vector opt.
369                  ! surface pressure gradient
370                  IF( lk_vvl) THEN
371                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
372                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
373                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
374                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
375                  ELSE
376                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
377                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
378                  ENDIF
379                  ! energy conserving formulation for planetary vorticity term
380                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
381                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
382                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
383                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
384                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
385                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
386                  ! after velocities with implicit bottom friction
387                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   &
388                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
389                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   &
390                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
391               END DO
392            END DO
393            !
394         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
395            DO jj = 2, jpjm1
396               DO ji = fs_2, fs_jpim1   ! vector opt.
397                   ! surface pressure gradient
398                  IF( lk_vvl) THEN
399                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
400                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
401                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
402                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
403                  ELSE
404                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
405                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
406                  ENDIF
407                  ! enstrophy conserving formulation for planetary vorticity term
408                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
409                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
410                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
411                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
412                  ! after velocities with implicit bottom friction
413                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   &
414                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
415                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   &
416                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
417               END DO
418            END DO
419            !
420         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
421            DO jj = 2, jpjm1
422               DO ji = fs_2, fs_jpim1   ! vector opt.
423                  ! surface pressure gradient
424                  IF( lk_vvl) THEN
425                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
426                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
427                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
428                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
429                  ELSE
430                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
431                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
432                  ENDIF
433                  ! energy/enstrophy conserving formulation for planetary vorticity term
434                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
435                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
436                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
437                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
438                  ! after velocities with implicit bottom friction
439                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   &
440                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
441                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   &
442                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
443               END DO
444            END DO
445            !
446         ENDIF
447         !                                                !* domain lateral boundary
448         !                                                !  -----------------------
449         !                                                      ! Flather's boundary condition for the barotropic loop :
450         !                                                      !         - Update sea surface height on each open boundary
451         !                                                      !         - Correct the velocity
452
453         IF( lk_obc                   )   CALL obc_fla_ts
454         IF( lk_bdy .OR. ln_bdy_tides )   CALL bdy_dyn_fla( sshn_e ) 
455         !
456         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
457         CALL lbc_lnk( va_e  , 'V', -1. )
458         CALL lbc_lnk( ssha_e, 'T',  1. )
459
460         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
461         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
462         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
463
464         !                                                !* Time filter and swap
465         !                                                !  --------------------
466         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
467            zsshb_e(:,:) = sshn_e(:,:)
468            zub_e  (:,:) = zun_e (:,:)
469            zvb_e  (:,:) = zvn_e (:,:)
470            sshn_e (:,:) = ssha_e(:,:)
471            zun_e  (:,:) = ua_e  (:,:)
472            zvn_e  (:,:) = va_e  (:,:)
473         ELSE                                                   ! Swap + Filter
474            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
475            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
476            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
477            sshn_e (:,:) = ssha_e(:,:)
478            zun_e  (:,:) = ua_e  (:,:)
479            zvn_e  (:,:) = va_e  (:,:)
480         ENDIF
481
482         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
483            !                                             !  ------------------
484            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
485               DO ji = 1, fs_jpim1   ! Vector opt.
486                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
487                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
488                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
489                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
490                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
491                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
492               END DO
493            END DO
494            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
495            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
496            !
497            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
498            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
499            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
500            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
501            !
502         ENDIF
503         !                                                 ! ==================== !
504      END DO                                               !        end loop      !
505      !                                                    ! ==================== !
506
507#if defined key_obc
508      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
509      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
510      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
511      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
512#endif
513
514      ! -----------------------------------------------------------------------------
515      ! Phase 3. update the general trend with the barotropic trend
516      ! -----------------------------------------------------------------------------
517      !
518      !                                   !* Time average ==> after barotropic u, v, ssh
519      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
520      un_b  (:,:) = zcoef * zu_sum  (:,:) 
521      vn_b  (:,:) = zcoef * zv_sum  (:,:) 
522      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
523      !
524      !                                   !* update the general momentum trend
525      DO jk=1,jpkm1
526         ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b
527         va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b
528      END DO
529      !
530      !                                   !* write time-spliting arrays in the restart
531      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
532      !
533      !
534   END SUBROUTINE dyn_spg_ts
535
536
537   SUBROUTINE ts_rst( kt, cdrw )
538      !!---------------------------------------------------------------------
539      !!                   ***  ROUTINE ts_rst  ***
540      !!
541      !! ** Purpose : Read or write time-splitting arrays in restart file
542      !!----------------------------------------------------------------------
543      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
544      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
545      !
546      INTEGER ::  ji, jk        ! dummy loop indices
547      !!----------------------------------------------------------------------
548      !
549      IF( TRIM(cdrw) == 'READ' ) THEN
550         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
551            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
552            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
553         ELSE
554            un_b (:,:) = 0.e0
555            vn_b (:,:) = 0.e0
556            ! vertical sum
557            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
558               DO jk = 1, jpkm1
559                  DO ji = 1, jpij
560                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
561                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
562                  END DO
563               END DO
564            ELSE                             ! No  vector opt.
565               DO jk = 1, jpkm1
566                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
567                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
568               END DO
569            ENDIF
570            un_b (:,:) = un_b(:,:) * hur(:,:)
571            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
572         ENDIF
573      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
574         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! external velocity issued
575         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
576      ENDIF
577      !
578   END SUBROUTINE ts_rst
579
580#else
581   !!----------------------------------------------------------------------
582   !!   Default case :   Empty module   No standart free surface cst volume
583   !!----------------------------------------------------------------------
584CONTAINS
585   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
586      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
587   END SUBROUTINE dyn_spg_ts
588   SUBROUTINE ts_rst( kt, cdrw )     ! Empty routine
589      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
590      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
591      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
592   END SUBROUTINE ts_rst   
593#endif
594   
595   !!======================================================================
596END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.