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

Last change on this file since 2724 was 2724, checked in by rblod, 10 years ago

Correct dynspg_ts, see ticket #812

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