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

Last change on this file since 2590 was 2590, checked in by trackstand2, 13 years ago

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