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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2382

Last change on this file since 2382 was 2338, checked in by rblod, 14 years ago

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