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/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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