New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynspg_ts.F90 in branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3057

Last change on this file since 3057 was 3057, checked in by hliu, 12 years ago

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