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/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2528

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

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