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

Last change on this file since 1779 was 1779, checked in by rblod, 12 years ago

Correct vector optimisation in dynspg_ts, see ticket #621

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