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

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