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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

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