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

Last change on this file since 2564 was 2564, checked in by rfurner, 13 years ago

bug fix for time splitting module

  • 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   USE zdf_oce
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC dyn_spg_ts  ! routine called by step.F90
48   PUBLIC ts_rst      ! routine called by istate.F90
49
50
51   REAL(wp), DIMENSION(jpi,jpj) ::  ftnw, ftne   ! triad of coriolis parameter
52   REAL(wp), DIMENSION(jpi,jpj) ::  ftsw, ftse   ! (only used with een vorticity scheme)
53
54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   un_b, vn_b   ! now    averaged velocity
55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ub_b, vb_b   ! before averaged velocity
56
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!-------------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
65   !!-------------------------------------------------------------------------
66
67CONTAINS
68
69   SUBROUTINE dyn_spg_ts( kt )
70      !!----------------------------------------------------------------------
71      !!                  ***  routine dyn_spg_ts  ***
72      !!
73      !! ** Purpose :   Compute the now trend due to the surface pressure
74      !!      gradient in case of free surface formulation with time-splitting.
75      !!      Add it to the general trend of momentum equation.
76      !!
77      !! ** Method  :   Free surface formulation with time-splitting
78      !!      -1- Save the vertically integrated trend. This general trend is
79      !!          held constant over the barotropic integration.
80      !!          The Coriolis force is removed from the general trend as the
81      !!          surface gradient and the Coriolis force are updated within
82      !!          the barotropic integration.
83      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and
84      !!          barotropic velocity (ua_e and va_e) through barotropic
85      !!          momentum and continuity integration. Barotropic former
86      !!          variables are time averaging over the full barotropic cycle
87      !!          (= 2 * baroclinic time step) and saved in uX_b
88      !!          and vX_b (X specifying after, now or before).
89      !!      -3- The new general trend becomes :
90      !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b )
91      !!
92      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
93      !!
94      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
95      !!---------------------------------------------------------------------
96      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
97      !!
98      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
99      INTEGER  ::   icycle           ! 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            ! Apply stability criteria for bottom friction
268            !RBbug for vvl and external mode we may need to use varying fse3
269            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
270            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
271            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
272         END DO
273      END DO
274
275      IF( lk_vvl ) THEN
276         DO jj = 2, jpjm1
277            DO ji = fs_2, fs_jpim1   ! vector opt.
278               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
279                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
280               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
281                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
282            END DO
283         END DO
284      ELSE
285         DO jj = 2, jpjm1
286            DO ji = fs_2, fs_jpim1   ! vector opt.
287               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
288               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
289            END DO
290         END DO
291      ENDIF
292
293      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
294      !                                   ! --------------------------
295      zua(:,:) = zua(:,:) * hur(:,:)
296      zva(:,:) = zva(:,:) * hvr(:,:)
297      !
298
299      ! -----------------------------------------------------------------------
300      !  Phase 2 : Integration of the barotropic equations with time splitting
301      ! -----------------------------------------------------------------------
302      !
303      !                                             ! ==================== !
304      !                                             !    Initialisations   !
305      !                                             ! ==================== !
306      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
307     
308      !                                ! Start from NOW field
309      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
310      hv_e   (:,:) = hv   (:,:) 
311      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
312      hvr_e  (:,:) = hvr  (:,:)
313!RBbug     zsshb_e(:,:) = sshn (:,:) 
314      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
315      sshn_e (:,:) = sshn (:,:)
316     
317      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
318      zvn_e  (:,:) = vn_b (:,:)
319      zub_e  (:,:) = un_b (:,:)
320      zvb_e  (:,:) = vn_b (:,:)
321
322      zu_sum  (:,:) = un_b (:,:)           ! summation
323      zv_sum  (:,:) = vn_b (:,:)
324      zssh_sum(:,:) = sshn (:,:)
325
326#if defined key_obc
327      ! set ssh corrections to 0
328      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
329      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
330      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
331      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
332      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
333#endif
334
335      !                                             ! ==================== !
336      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
337         !                                          ! ==================== !
338         z2dt_e = 2. * ( rdt / nn_baro )
339         IF( jn == 1 )   z2dt_e = rdt / nn_baro
340
341         !                                                !* Update the forcing (OBC, BDY and tides)
342         !                                                !  ------------------
343         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
344         IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle )
345
346         !                                                !* after ssh_e
347         !                                                !  -----------
348         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
349            DO ji = fs_2, fs_jpim1   ! vector opt.
350               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
351                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
352                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
353                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
354            END DO
355         END DO
356         !
357#if defined key_obc
358         !                                                     ! OBC : zhdiv must be zero behind the open boundary
359!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
360         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
361         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
362         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
363         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
364#endif
365#if defined key_bdy
366         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
367#endif
368         !
369         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
370            DO ji = fs_2, fs_jpim1   ! vector opt.
371               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
372            END DO
373         END DO
374
375         !                                                !* after barotropic velocities (vorticity scheme dependent)
376         !                                                !  --------------------------- 
377         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
378         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
379         !
380         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
381            DO jj = 2, jpjm1
382               DO ji = fs_2, fs_jpim1   ! vector opt.
383                  ! surface pressure gradient
384                  IF( lk_vvl) THEN
385                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
386                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
387                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
388                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
389                  ELSE
390                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
391                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
392                  ENDIF
393                  ! energy conserving formulation for planetary vorticity term
394                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
395                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
396                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
397                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
398                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
399                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
400                  ! after velocities with implicit bottom friction
401                  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)   &
402                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
403                  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)   &
404                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
405               END DO
406            END DO
407            !
408         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
409            DO jj = 2, jpjm1
410               DO ji = fs_2, fs_jpim1   ! vector opt.
411                   ! surface pressure gradient
412                  IF( lk_vvl) THEN
413                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
414                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
415                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
416                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
417                  ELSE
418                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
419                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
420                  ENDIF
421                  ! enstrophy conserving formulation for planetary vorticity term
422                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
423                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
424                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
425                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
426                  ! after velocities with implicit bottom friction
427                  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)   &
428                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
429                  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)   &
430                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
431               END DO
432            END DO
433            !
434         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
435            DO jj = 2, jpjm1
436               DO ji = fs_2, fs_jpim1   ! vector opt.
437                  ! surface pressure gradient
438                  IF( lk_vvl) THEN
439                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
440                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
441                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
442                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
443                  ELSE
444                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
445                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
446                  ENDIF
447                  ! energy/enstrophy conserving formulation for planetary vorticity term
448                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
449                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
450                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
451                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
452                  ! after velocities with implicit bottom friction
453                  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)   &
454                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
455                  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)   &
456                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
457               END DO
458            END DO
459            !
460         ENDIF
461         !                                                !* domain lateral boundary
462         !                                                !  -----------------------
463         !                                                      ! Flather's boundary condition for the barotropic loop :
464         !                                                      !         - Update sea surface height on each open boundary
465         !                                                      !         - Correct the velocity
466
467         IF( lk_obc               )   CALL obc_fla_ts
468         IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e ) 
469         !
470         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
471         CALL lbc_lnk( va_e  , 'V', -1. )
472         CALL lbc_lnk( ssha_e, 'T',  1. )
473
474         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
475         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
476         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
477
478         !                                                !* Time filter and swap
479         !                                                !  --------------------
480         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
481            zsshb_e(:,:) = sshn_e(:,:)
482            zub_e  (:,:) = zun_e (:,:)
483            zvb_e  (:,:) = zvn_e (:,:)
484            sshn_e (:,:) = ssha_e(:,:)
485            zun_e  (:,:) = ua_e  (:,:)
486            zvn_e  (:,:) = va_e  (:,:)
487         ELSE                                                   ! Swap + Filter
488            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
489            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
490            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
491            sshn_e (:,:) = ssha_e(:,:)
492            zun_e  (:,:) = ua_e  (:,:)
493            zvn_e  (:,:) = va_e  (:,:)
494         ENDIF
495
496         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
497            !                                             !  ------------------
498            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
499               DO ji = 1, fs_jpim1   ! Vector opt.
500                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
501                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
502                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
503                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
504                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
505                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
506               END DO
507            END DO
508            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
509            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
510            !
511            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
512            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
513            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
514            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
515            !
516         ENDIF
517         !                                                 ! ==================== !
518      END DO                                               !        end loop      !
519      !                                                    ! ==================== !
520
521#if defined key_obc
522      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
523      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
524      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
525      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
526#endif
527
528      ! -----------------------------------------------------------------------------
529      ! Phase 3. update the general trend with the barotropic trend
530      ! -----------------------------------------------------------------------------
531      !
532      !                                   !* Time average ==> after barotropic u, v, ssh
533      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
534      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
535      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
536      !
537      !                                   !* update the general momentum trend
538      DO jk=1,jpkm1
539         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
540         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
541      END DO
542      ub_b (:,:) = un_b(:,:)
543      vb_b (:,:) = vn_b(:,:)
544      un_b  (:,:) =  zu_sum(:,:) 
545      vn_b  (:,:) =  zv_sum(:,:) 
546      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
547      !
548      !                                   !* write time-spliting arrays in the restart
549      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
550      !
551      !
552   END SUBROUTINE dyn_spg_ts
553
554
555   SUBROUTINE ts_rst( kt, cdrw )
556      !!---------------------------------------------------------------------
557      !!                   ***  ROUTINE ts_rst  ***
558      !!
559      !! ** Purpose : Read or write time-splitting arrays in restart file
560      !!----------------------------------------------------------------------
561      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
562      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
563      !
564      INTEGER ::  ji, jk        ! dummy loop indices
565      !!----------------------------------------------------------------------
566      !
567      IF( TRIM(cdrw) == 'READ' ) THEN
568         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
569            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
570            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
571         ELSE
572            un_b (:,:) = 0.e0
573            vn_b (:,:) = 0.e0
574            ! vertical sum
575            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
576               DO jk = 1, jpkm1
577                  DO ji = 1, jpij
578                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
579                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
580                  END DO
581               END DO
582            ELSE                             ! No  vector opt.
583               DO jk = 1, jpkm1
584                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
585                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
586               END DO
587            ENDIF
588            un_b (:,:) = un_b(:,:) * hur(:,:)
589            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
590         ENDIF
591
592         ! Vertically integrated velocity (before)
593         IF (neuler/=0) THEN
594            ub_b (:,:) = 0.e0
595            vb_b (:,:) = 0.e0
596
597            ! vertical sum
598            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
599               DO jk = 1, jpkm1
600                  DO ji = 1, jpij
601                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
602                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
603                  END DO
604               END DO
605            ELSE                             ! No  vector opt.
606               DO jk = 1, jpkm1
607                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
608                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
609               END DO
610            ENDIF
611
612            IF( lk_vvl ) THEN
613               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
614               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
615            ELSE
616               ub_b(:,:) = ub_b(:,:) * hur(:,:)
617               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
618            ENDIF
619         ELSE                                 ! neuler==0
620            ub_b (:,:) = un_b (:,:)
621            vb_b (:,:) = vn_b (:,:)
622         ENDIF
623
624         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
625            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered extrenal ssh
626         ELSE
627            sshn_b(:,:)=sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
628         ENDIF
629      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
630         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
631         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
632         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !
633      ENDIF
634      !
635   END SUBROUTINE ts_rst
636
637#else
638   !!----------------------------------------------------------------------
639   !!   Default case :   Empty module   No standart free surface cst volume
640   !!----------------------------------------------------------------------
641CONTAINS
642   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
643      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
644   END SUBROUTINE dyn_spg_ts
645   SUBROUTINE ts_rst( kt, cdrw )     ! Empty routine
646      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
647      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
648      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
649   END SUBROUTINE ts_rst   
650#endif
651   
652   !!======================================================================
653END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.