New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2865

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