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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2684

Last change on this file since 2684 was 2684, checked in by rfurner, 13 years ago

bug fixes for running with key_bdy
  • Property svn:keywords set to Id
File size: 36.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  !!---------------------------------------------------------------------
11#if defined key_dynspg_ts   ||   defined key_esopa
12   !!----------------------------------------------------------------------
13   !!   'key_dynspg_ts'         free surface cst volume with time splitting
14   !!----------------------------------------------------------------------
15   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
16   !!                 splitting scheme and add to the general trend
17   !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! surface boundary condition: ocean
22   USE dynspg_oce      ! surface pressure gradient variables
23   USE phycst          ! physical constants
24   USE domvvl          ! variable volume
25   USE zdfbfr          ! bottom friction
26   USE obcdta          ! open boundary condition data     
27   USE obcfla          ! Flather open boundary condition 
28   USE dynvor          ! vorticity term
29   USE obc_oce         ! Lateral open boundary condition
30   USE obc_par         ! open boundary condition parameters
31   USE bdy_oce         ! unstructured open boundaries
32   USE bdy_par         ! unstructured open boundaries
33   USE bdydta          ! unstructured open boundaries
34   USE bdydyn          ! unstructured open boundaries
35   USE bdytides        ! tidal forcing at unstructured open boundaries.
36   USE lib_mpp         ! distributed memory computing library
37   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
38   USE prtctl          ! Print control
39   USE in_out_manager  ! I/O manager
40   USE iom             ! IOM library
41   USE restart         ! only for lrst_oce
42   USE zdf_oce
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC dyn_spg_ts        ! routine called by step.F90
48   PUBLIC ts_rst            ! routine called by istate.F90
49   PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90
50
51
52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter
53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme)
54
55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b       ! now    averaged velocity
56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b       ! before averaged velocity
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   INTEGER FUNCTION dyn_spg_ts_alloc()
69      !!----------------------------------------------------------------------
70      !!                  ***  routine dyn_spg_ts_alloc  ***
71      !!----------------------------------------------------------------------
72      !
73      ALLOCATE( ftnw  (jpi,jpj) , ftne  (jpi,jpj) , ftsw  (jpi,jpj) , ftse (jpi,jpj) ,                                       &
74         &      un_b  (jpi,jpj) , vn_b  (jpi,jpj) , ub_b  (jpi,jpj) , vb_b (jpi,jpj) , ua_e  (jpi,jpj) , va_e  (jpi,jpj) ,   &
75         &      sshn_e(jpi,jpj) , ssha_e(jpi,jpj) , sshn_b(jpi,jpj) ,                                                        &
76         &      hu_e  (jpi,jpj) , hv_e  (jpi,jpj) , hur_e (jpi,jpj) , hvr_e(jpi,jpj) , STAT=dyn_spg_ts_alloc )
77      IF(lk_mpp)   CALL mpp_sum( dyn_spg_ts_alloc )
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: zhdiv => wrk_2d_1, zsshb_e => wrk_2d_2
111      USE wrk_nemo, ONLY: zbfru => wrk_2d_3  , zbfrv => wrk_2d_4 
112      USE wrk_nemo, ONLY: zsshun_e => wrk_2d_5, zsshvn_e => wrk_2d_6
113      USE wrk_nemo, ONLY: zcu => wrk_2d_7, zwx => wrk_2d_8, zua => wrk_2d_9, zun => wrk_2d_10
114      USE wrk_nemo, ONLY: zcv => wrk_2d_11, zwy => wrk_2d_12, zva => wrk_2d_13, zvn => wrk_2d_14
115      USE wrk_nemo, ONLY: zun_e => wrk_2d_15, zub_e => wrk_2d_16, zu_sum => wrk_2d_17
116      USE wrk_nemo, ONLY: zvn_e => wrk_2d_18, zvb_e => wrk_2d_19, zv_sum => wrk_2d_20
117      USE wrk_nemo, ONLY: zssh_sum => wrk_2d_21
118      !!
119      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
120      !!
121      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
122      INTEGER  ::   icycle           ! local scalar
123      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars
124      REAL(wp) ::   z1_8, zx1, zy1                   !   -      -
125      REAL(wp) ::   z1_4, zx2, zy2                   !   -      -
126      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !   -      -
127      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !   -      -
128      !!----------------------------------------------------------------------
129
130      IF( wrk_in_use(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,        &
131         &              11,12,13,14,15,16,17,18,19,20,21 ) ) THEN
132         CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' )   ;   RETURN
133      END IF
134
135      IF( kt == nit000 ) THEN             !* initialisation
136         !
137         IF(lwp) WRITE(numout,*)
138         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
139         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
140         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
141         !
142         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b 
143         !
144         ua_e  (:,:) = un_b (:,:)
145         va_e  (:,:) = vn_b (:,:)
146         hu_e  (:,:) = hu   (:,:)
147         hv_e  (:,:) = hv   (:,:)
148         hur_e (:,:) = hur  (:,:)
149         hvr_e (:,:) = hvr  (:,:)
150         IF( ln_dynvor_een ) THEN
151            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
152            DO jj = 2, jpj
153               DO ji = fs_2, jpi   ! vector opt.
154                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
155                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
156                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
157                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
158               END DO
159            END DO
160         ENDIF
161         !
162      ENDIF
163
164      !                                   !* Local constant initialization
165      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
166      z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
167      z1_4 = 0.5 * 0.5
168      zraur  = 1. / rau0                                    ! 1 / volumic mass
169      !
170      zhdiv(:,:) = 0.e0                                     ! barotropic divergence
171      zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F)
172      zv_sld = 0.e0   ;   zv_asp = 0.e0
173
174      ! -----------------------------------------------------------------------------
175      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
176      ! -----------------------------------------------------------------------------
177      !     
178      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
179      !                                   ! --------------------------
180      zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0 
181      zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0 
182      !
183      DO jk = 1, jpkm1
184#if defined key_vectopt_loop
185         DO jj = 1, 1         !Vector opt. => forced unrolling
186            DO ji = 1, jpij
187#else
188         DO jj = 1, jpj
189            DO ji = 1, jpi
190#endif
191               !                                                                              ! now trend
192               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
193               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
194               !                                                                              ! now velocity
195               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
196               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
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
314      ! -----------------------------------------------------------------------
315      !  Phase 2 : Integration of the barotropic equations with time splitting
316      ! -----------------------------------------------------------------------
317      !
318      !                                             ! ==================== !
319      !                                             !    Initialisations   !
320      !                                             ! ==================== !
321      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
322     
323      !                                ! Start from NOW field
324      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
325      hv_e   (:,:) = hv   (:,:) 
326      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
327      hvr_e  (:,:) = hvr  (:,:)
328!RBbug     zsshb_e(:,:) = sshn (:,:) 
329      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
330      sshn_e (:,:) = sshn (:,:)
331     
332      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
333      zvn_e  (:,:) = vn_b (:,:)
334      zub_e  (:,:) = un_b (:,:)
335      zvb_e  (:,:) = vn_b (:,:)
336
337      zu_sum  (:,:) = un_b (:,:)           ! summation
338      zv_sum  (:,:) = vn_b (:,:)
339      zssh_sum(:,:) = sshn (:,:)
340
341#if defined key_obc
342      ! set ssh corrections to 0
343      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
344      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
345      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
346      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
347      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
348#endif
349
350      !                                             ! ==================== !
351      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
352         !                                          ! ==================== !
353         z2dt_e = 2. * ( rdt / nn_baro )
354         IF( jn == 1 )   z2dt_e = rdt / nn_baro
355
356         !                                                !* Update the forcing (OBC, BDY and tides)
357         !                                                !  ------------------
358         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
359         IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle )
360
361         !                                                !* after ssh_e
362         !                                                !  -----------
363         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
364            DO ji = fs_2, fs_jpim1   ! vector opt.
365               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
366                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
367                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
368                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
369            END DO
370         END DO
371         !
372#if defined key_obc
373         !                                                     ! OBC : zhdiv must be zero behind the open boundary
374!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
375         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
376         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
377         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
378         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
379#endif
380#if defined key_bdy
381         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
382#endif
383         !
384         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
385            DO ji = fs_2, fs_jpim1   ! vector opt.
386               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
387            END DO
388         END DO
389
390         !                                                !* after barotropic velocities (vorticity scheme dependent)
391         !                                                !  --------------------------- 
392         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
393         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
394         !
395         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
396            DO jj = 2, jpjm1
397               DO ji = fs_2, fs_jpim1   ! vector opt.
398                  ! surface pressure gradient
399                  IF( lk_vvl) THEN
400                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
401                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
402                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
403                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
404                  ELSE
405                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
406                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
407                  ENDIF
408                  ! energy conserving formulation for planetary vorticity term
409                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
410                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
411                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
412                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
413                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
414                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
415                  ! after velocities with implicit bottom friction
416                  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)   &
417                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
418                  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)   &
419                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
420               END DO
421            END DO
422            !
423         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
424            DO jj = 2, jpjm1
425               DO ji = fs_2, fs_jpim1   ! vector opt.
426                   ! surface pressure gradient
427                  IF( lk_vvl) THEN
428                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    &
429                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
430                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
431                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
432                  ELSE
433                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
434                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
435                  ENDIF
436                  ! enstrophy conserving formulation for planetary vorticity term
437                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
438                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
439                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
440                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
441                  ! after velocities with implicit bottom friction
442                  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)   &
443                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
444                  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)   &
445                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
446               END DO
447            END DO
448            !
449         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
450            DO jj = 2, jpjm1
451               DO ji = fs_2, fs_jpim1   ! vector opt.
452                  ! surface pressure gradient
453                  IF( lk_vvl) THEN
454                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    &
455                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj)
456                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    &
457                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj)
458                  ELSE
459                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
460                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
461                  ENDIF
462                  ! energy/enstrophy conserving formulation for planetary vorticity term
463                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
464                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
465                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
466                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
467                  ! after velocities with implicit bottom friction
468                  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)   &
469                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
470                  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)   &
471                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
472               END DO
473            END DO
474            !
475         ENDIF
476         !                                                !* domain lateral boundary
477         !                                                !  -----------------------
478         !                                                      ! Flather's boundary condition for the barotropic loop :
479         !                                                      !         - Update sea surface height on each open boundary
480         !                                                      !         - Correct the velocity
481
482         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
483         IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e ) 
484         !
485         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
486         CALL lbc_lnk( va_e  , 'V', -1. )
487         CALL lbc_lnk( ssha_e, 'T',  1. )
488
489         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
490         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
491         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
492
493         !                                                !* Time filter and swap
494         !                                                !  --------------------
495         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
496            zsshb_e(:,:) = sshn_e(:,:)
497            zub_e  (:,:) = zun_e (:,:)
498            zvb_e  (:,:) = zvn_e (:,:)
499            sshn_e (:,:) = ssha_e(:,:)
500            zun_e  (:,:) = ua_e  (:,:)
501            zvn_e  (:,:) = va_e  (:,:)
502         ELSE                                                   ! Swap + Filter
503            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
504            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
505            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
506            sshn_e (:,:) = ssha_e(:,:)
507            zun_e  (:,:) = ua_e  (:,:)
508            zvn_e  (:,:) = va_e  (:,:)
509         ENDIF
510
511         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
512            !                                             !  ------------------
513            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
514               DO ji = 1, fs_jpim1   ! Vector opt.
515                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
516                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
517                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
518                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
519                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
520                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
521               END DO
522            END DO
523            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
524            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
525            !
526            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
527            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
528            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
529            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
530            !
531         ENDIF
532         !                                                 ! ==================== !
533      END DO                                               !        end loop      !
534      !                                                    ! ==================== !
535
536#if defined key_obc
537      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
538      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
539      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
540      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
541#endif
542
543      ! -----------------------------------------------------------------------------
544      ! Phase 3. update the general trend with the barotropic trend
545      ! -----------------------------------------------------------------------------
546      !
547      !                                   !* Time average ==> after barotropic u, v, ssh
548      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
549      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
550      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
551      !
552      !                                   !* update the general momentum trend
553      DO jk=1,jpkm1
554         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
555         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
556      END DO
557      ub_b (:,:) = un_b(:,:)
558      vb_b (:,:) = vn_b(:,:)
559      un_b  (:,:) =  zu_sum(:,:) 
560      vn_b  (:,:) =  zv_sum(:,:) 
561      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
562      !
563      !                                   !* write time-spliting arrays in the restart
564      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
565      !
566      !
567      IF( wrk_not_released(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,         &
568         &                    11,12,13,14,15,16,17,18,19,20,21) )    &
569         CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays.')
570      !
571   END SUBROUTINE dyn_spg_ts
572
573
574   SUBROUTINE ts_rst( kt, cdrw )
575      !!---------------------------------------------------------------------
576      !!                   ***  ROUTINE ts_rst  ***
577      !!
578      !! ** Purpose : Read or write time-splitting arrays in restart file
579      !!----------------------------------------------------------------------
580      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
581      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
582      !
583      INTEGER ::  ji, jk        ! dummy loop indices
584      !!----------------------------------------------------------------------
585      !
586      IF( TRIM(cdrw) == 'READ' ) THEN
587         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
588            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
589            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
590         ELSE
591            un_b (:,:) = 0.e0
592            vn_b (:,:) = 0.e0
593            ! vertical sum
594            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
595               DO jk = 1, jpkm1
596                  DO ji = 1, jpij
597                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
598                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
599                  END DO
600               END DO
601            ELSE                             ! No  vector opt.
602               DO jk = 1, jpkm1
603                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
604                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
605               END DO
606            ENDIF
607            un_b (:,:) = un_b(:,:) * hur(:,:)
608            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
609         ENDIF
610
611         ! Vertically integrated velocity (before)
612         IF (neuler/=0) THEN
613            ub_b (:,:) = 0.e0
614            vb_b (:,:) = 0.e0
615
616            ! vertical sum
617            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
618               DO jk = 1, jpkm1
619                  DO ji = 1, jpij
620                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
621                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
622                  END DO
623               END DO
624            ELSE                             ! No  vector opt.
625               DO jk = 1, jpkm1
626                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
627                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
628               END DO
629            ENDIF
630
631            IF( lk_vvl ) THEN
632               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
633               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
634            ELSE
635               ub_b(:,:) = ub_b(:,:) * hur(:,:)
636               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
637            ENDIF
638         ELSE                                 ! neuler==0
639            ub_b (:,:) = un_b (:,:)
640            vb_b (:,:) = vn_b (:,:)
641         ENDIF
642
643         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
644            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered extrenal ssh
645         ELSE
646            sshn_b(:,:)=sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
647         ENDIF
648      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
649         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
650         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
651         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !
652      ENDIF
653      !
654   END SUBROUTINE ts_rst
655
656#else
657   !!----------------------------------------------------------------------
658   !!   Default case :   Empty module   No standart free surface cst volume
659   !!----------------------------------------------------------------------
660CONTAINS
661   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
662      dyn_spg_ts_alloc = 0
663   END FUNCTION dyn_spg_ts_alloc
664   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
665      INTEGER, INTENT(in) :: kt
666      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
667   END SUBROUTINE dyn_spg_ts
668   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
669      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
670      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
671      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
672   END SUBROUTINE ts_rst   
673#endif
674   
675   !!======================================================================
676END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.