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

source: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3008

Last change on this file since 3008 was 3008, checked in by acc, 13 years ago

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