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

Last change on this file since 2872 was 2872, checked in by hliu, 13 years ago

1) added the semi-implicit bottom friction code (3 in DYN, 1 in ZDF), 2) added par_AMM7/12.h90, 3) added two arch files for NOCL mobius and ubuntu linux
  • Property svn:keywords set to Id
File size: 38.5 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     ! 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      ! -----------------------------------------------------------------------------
175      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
176      ! -----------------------------------------------------------------------------
177      !     
178      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
179      !                                   ! --------------------------
180      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp
181      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
182      !
183      DO jk = 1, jpkm1
184#if defined key_vectopt_loop
185         DO jj = 1, 1         !Vector opt. => forced unrolling
186            DO ji = 1, jpij
187#else
188         DO jj = 1, jpj
189            DO ji = 1, jpi
190#endif
191               !                                                                              ! now trend
192               zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
193               zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
194               !                                                                              ! now velocity
195               zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk)
196               zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)               
197               !
198#if defined key_vvl
199               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) 
200               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)   
201#else
202               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
203               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
204#endif
205            END DO
206         END DO
207      END DO
208
209
210      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent)
211      !                                   ! ---------------------------====
212      zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport
213      zwy(:,:) = zvn(:,:) * e1v(:,:)
214      !
215      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
219               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
220               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
221               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
222               ! energy conserving formulation for planetary vorticity term
223               zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
224               zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
225            END DO
226         END DO
227         !
228      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
229         DO jj = 2, jpjm1
230            DO ji = fs_2, fs_jpim1   ! vector opt.
231               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
232               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
233               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
234               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
235            END DO
236         END DO
237         !
238      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
239         DO jj = 2, jpjm1
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
242                  &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
243               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)   &
244                  &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
245            END DO
246         END DO
247         !
248      ENDIF
249
250      !                                   !* Right-Hand-Side of the barotropic momentum equation
251      !                                   ! ----------------------------------------------------
252      IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient
253         DO jj = 2, jpjm1 
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255      !     remove the rhd term according to J. Chanut's suggestion           
256               zcu(ji,jj) = zcu(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj)
257               zcv(ji,jj) = zcv(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj)
258            END DO
259         END DO
260      ENDIF
261
262     IF(ln_bfrimp.AND.lk_vvl) THEN
263!    Remove the bottom stress trend from 3-D sea surface level gradient and Coriolis force           
264!    H. Liu, Sept. 2011
265      DO jj = 2, jpjm1         
266         DO ji = fs_2, fs_jpim1
267            ikbu = mbku(ji,jj)
268            ikbv = mbkv(ji,jj)
269            ua_btm = zcu(ji,jj) * z2dt_b * hur(ji,jj) * umask (ji,jj,ikbu)
270            va_btm = zcv(ji,jj) * z2dt_b * hvr(ji,jj) * vmask (ji,jj,ikbv)
271
272            ua(ji,jj,ikbu) = ua(ji,jj,ikbu) - bfrua(ji,jj) * ua_btm / fse3u(ji,jj,ikbu)
273            va(ji,jj,ikbv) = va(ji,jj,ikbv) - bfrva(ji,jj) * va_btm / fse3v(ji,jj,ikbv)
274           
275            zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
276            zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
277         END DO
278      END DO
279     END IF
280      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
281      DO jk = 1, jpkm1                    ! --------------------------
282         DO jj = 2, jpjm1
283            DO ji = fs_2, fs_jpim1   ! vector opt.
284               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
285               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
286            END DO
287         END DO
288      END DO
289
290      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend
291         DO ji = fs_2, fs_jpim1
292            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
293            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
294         END DO
295      END DO
296
297                   
298      !                                             ! Remove barotropic contribution of bottom friction
299      !                                             ! from the barotropic transport trend
300      zcoef = -1. / z2dt_b
301      IF(.NOT.ln_bfrimp) THEN
302# if defined key_vectopt_loop
303      DO jj = 1, 1
304         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
305# else
306      DO jj = 2, jpjm1
307         DO ji = 2, jpim1
308# endif
309            ! Apply stability criteria for bottom friction
310            !RBbug for vvl and external mode we may need to use varying fse3
311            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
312            zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  )
313            zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  )
314         END DO
315      END DO
316
317      IF( lk_vvl ) THEN
318         DO jj = 2, jpjm1
319            DO ji = fs_2, fs_jpim1   ! vector opt.
320               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   &
321                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
322               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   &
323                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
324            END DO
325         END DO
326      ELSE
327         DO jj = 2, jpjm1
328            DO ji = fs_2, fs_jpim1   ! vector opt.
329               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
330               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
331            END DO
332         END DO
333      ENDIF
334      END IF
335
336      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
337      !                                   ! --------------------------
338      zua(:,:) = zua(:,:) * hur(:,:)
339      zva(:,:) = zva(:,:) * hvr(:,:)
340      !
341      IF( lk_vvl ) THEN
342         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
343         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
344      ELSE
345         ub_b(:,:) = ub_b(:,:) * hur(:,:)
346         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
347      ENDIF
348
349      ! -----------------------------------------------------------------------
350      !  Phase 2 : Integration of the barotropic equations with time splitting
351      ! -----------------------------------------------------------------------
352      !
353      !                                             ! ==================== !
354      !                                             !    Initialisations   !
355      !                                             ! ==================== !
356      icycle = 2  * nn_baro            ! Number of barotropic sub time-step
357     
358      !                                ! Start from NOW field
359      hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points
360      hv_e   (:,:) = hv   (:,:) 
361      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points
362      hvr_e  (:,:) = hvr  (:,:)
363!RBbug     zsshb_e(:,:) = sshn (:,:) 
364      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now)
365      sshn_e (:,:) = sshn (:,:)
366     
367      zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external)
368      zvn_e  (:,:) = vn_b (:,:)
369      zub_e  (:,:) = un_b (:,:)
370      zvb_e  (:,:) = vn_b (:,:)
371
372      zu_sum  (:,:) = un_b (:,:)           ! summation
373      zv_sum  (:,:) = vn_b (:,:)
374      zssh_sum(:,:) = sshn (:,:)
375
376#if defined key_obc
377      ! set ssh corrections to 0
378      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
379      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
380      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
381      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
382      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
383#endif
384
385      !                                             ! ==================== !
386      DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1)
387         !                                          ! ==================== !
388         z2dt_e = 2. * ( rdt / nn_baro )
389         IF( jn == 1 )   z2dt_e = rdt / nn_baro
390
391         !                                                !* Update the forcing (OBC, BDY and tides)
392         !                                                !  ------------------
393         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
394         IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle )
395
396         !                                                !* after ssh_e
397         !                                                !  -----------
398         DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports
399            DO ji = fs_2, fs_jpim1   ! vector opt.
400               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     &
401                  &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     &
402                  &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     &
403                  &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) )
404            END DO
405         END DO
406         !
407#if defined key_obc
408         !                                                     ! OBC : zhdiv must be zero behind the open boundary
409!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
410         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
411         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
412         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
413         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
414#endif
415#if defined key_bdy
416         zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask
417#endif
418         !
419         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e
420            DO ji = fs_2, fs_jpim1   ! vector opt.
421               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 
422            END DO
423         END DO
424
425         !                                                !* after barotropic velocities (vorticity scheme dependent)
426         !                                                !  --------------------------- 
427         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport
428         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
429         !
430         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==!
431            DO jj = 2, jpjm1
432               DO ji = fs_2, fs_jpim1   ! vector opt.
433                  ! surface pressure gradient
434                  zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
435                  zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
436                  ! energy conserving formulation for planetary vorticity term
437                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
438                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
439                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
440                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
441                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
442                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
443                  ! after velocities with implicit bottom friction
444
445                  IF(ln_bfrimp) THEN
446                  !   A new method to implement the implicit bottom friction.
447                  !   H. Liu
448                  !   May 2010
449                  ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   &
450                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
451                  ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
452                  va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   &
453                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
454                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
455                  ELSE
456                  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)   &
457                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
458                  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)   &
459                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
460                  ENDIF
461               END DO
462            END DO
463            !
464         ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==!
465            DO jj = 2, jpjm1
466               DO ji = fs_2, fs_jpim1   ! vector opt.
467                   ! surface pressure gradient
468                  zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
469                  zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
470                  ! enstrophy conserving formulation for planetary vorticity term
471                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
472                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
473                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
474                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj)
475                  ! after velocities with implicit bottom friction
476                  IF(ln_bfrimp) THEN
477                  !   A new method to implement the implicit bottom friction.
478                  !   H. Liu
479                  !   May 2010
480                  ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   &
481                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
482                  ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
483                  va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   &
484                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
485                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
486                  ELSE
487                  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)   &
488                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
489                  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)   &
490                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
491                  ENDIF
492               END DO
493            END DO
494            !
495         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
496            DO jj = 2, jpjm1
497               DO ji = fs_2, fs_jpim1   ! vector opt.
498                  ! surface pressure gradient
499                  zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
500                  zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
501                  ! energy/enstrophy conserving formulation for planetary vorticity term
502                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
503                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
504                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   &
505                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj)
506                  ! after velocities with implicit bottom friction
507                  IF(ln_bfrimp) THEN
508                  !   A new method to implement the implicit bottom friction.
509                  !   H. Liu
510                  !   May 2010
511                  ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   &
512                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
513                  ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)   
514                  va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   &
515                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
516                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)   
517                  ELSE
518                  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)   &
519                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
520                  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)   &
521                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
522                  ENDIF
523               END DO
524            END DO
525            !
526         ENDIF
527         !                                                !* domain lateral boundary
528         !                                                !  -----------------------
529         !                                                      ! Flather's boundary condition for the barotropic loop :
530         !                                                      !         - Update sea surface height on each open boundary
531         !                                                      !         - Correct the velocity
532
533         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
534         IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e ) 
535         !
536         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries
537         CALL lbc_lnk( va_e  , 'V', -1. )
538         CALL lbc_lnk( ssha_e, 'T',  1. )
539
540         zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps
541         zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:) 
542         zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 
543
544         !                                                !* Time filter and swap
545         !                                                !  --------------------
546         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step)
547            zsshb_e(:,:) = sshn_e(:,:)
548            zub_e  (:,:) = zun_e (:,:)
549            zvb_e  (:,:) = zvn_e (:,:)
550            sshn_e (:,:) = ssha_e(:,:)
551            zun_e  (:,:) = ua_e  (:,:)
552            zvn_e  (:,:) = va_e  (:,:)
553         ELSE                                                   ! Swap + Filter
554            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
555            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:)
556            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:)
557            sshn_e (:,:) = ssha_e(:,:)
558            zun_e  (:,:) = ua_e  (:,:)
559            zvn_e  (:,:) = va_e  (:,:)
560         ENDIF
561
562         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only)
563            !                                             !  ------------------
564            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points
565               DO ji = 1, fs_jpim1   ! Vector opt.
566                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       &
567                     &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
568                     &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
569                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       &
570                     &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
571                     &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
572               END DO
573            END DO
574            CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions
575            CALL lbc_lnk( zsshvn_e, 'V', 1. ) 
576            !
577            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points
578            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
579            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
580            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
581            !
582         ENDIF
583         !                                                 ! ==================== !
584      END DO                                               !        end loop      !
585      !                                                    ! ==================== !
586
587#if defined key_obc
588      IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ?????
589      IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
590      IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
591      IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
592#endif
593
594      ! -----------------------------------------------------------------------------
595      ! Phase 3. update the general trend with the barotropic trend
596      ! -----------------------------------------------------------------------------
597      !
598      !                                   !* Time average ==> after barotropic u, v, ssh
599      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
600      zu_sum(:,:) = zcoef * zu_sum  (:,:) 
601      zv_sum(:,:) = zcoef * zv_sum  (:,:) 
602      !
603      !                                   !* update the general momentum trend
604      DO jk=1,jpkm1
605         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
606         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
607      END DO
608      un_b  (:,:) =  zu_sum(:,:) 
609      vn_b  (:,:) =  zv_sum(:,:) 
610      sshn_b(:,:) = zcoef * zssh_sum(:,:) 
611      !
612      !                                   !* write time-spliting arrays in the restart
613      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
614      !
615      !
616      IF( wrk_not_released(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,         &
617         &                    11,12,13,14,15,16,17,18,19,20,21) )    &
618         CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays')
619      !
620   END SUBROUTINE dyn_spg_ts
621
622
623   SUBROUTINE ts_rst( kt, cdrw )
624      !!---------------------------------------------------------------------
625      !!                   ***  ROUTINE ts_rst  ***
626      !!
627      !! ** Purpose : Read or write time-splitting arrays in restart file
628      !!----------------------------------------------------------------------
629      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
630      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
631      !
632      INTEGER ::  ji, jk        ! dummy loop indices
633      !!----------------------------------------------------------------------
634      !
635      IF( TRIM(cdrw) == 'READ' ) THEN
636         IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
637            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued
638            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
639         ELSE
640            un_b (:,:) = 0.e0
641            vn_b (:,:) = 0.e0
642            ! vertical sum
643            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
644               DO jk = 1, jpkm1
645                  DO ji = 1, jpij
646                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
647                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
648                  END DO
649               END DO
650            ELSE                             ! No  vector opt.
651               DO jk = 1, jpkm1
652                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
653                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
654               END DO
655            ENDIF
656            un_b (:,:) = un_b(:,:) * hur(:,:)
657            vn_b (:,:) = vn_b(:,:) * hvr(:,:)
658         ENDIF
659
660         ! Vertically integrated velocity (before)
661         IF (neuler/=0) THEN
662            ub_b (:,:) = 0.e0
663            vb_b (:,:) = 0.e0
664
665            ! vertical sum
666            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
667               DO jk = 1, jpkm1
668                  DO ji = 1, jpij
669                     ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
670                     vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
671                  END DO
672               END DO
673            ELSE                             ! No  vector opt.
674               DO jk = 1, jpkm1
675                  ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
676                  vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
677               END DO
678            ENDIF
679
680            IF( lk_vvl ) THEN
681               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
682               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
683            ELSE
684               ub_b(:,:) = ub_b(:,:) * hur(:,:)
685               vb_b(:,:) = vb_b(:,:) * hvr(:,:)
686            ENDIF
687         ELSE                                 ! neuler==0
688            ub_b (:,:) = un_b (:,:)
689            vb_b (:,:) = vn_b (:,:)
690         ENDIF
691
692         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
693            CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh
694         ELSE
695            sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value   
696         ENDIF
697      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
698         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh
699         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop
700         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !
701      ENDIF
702      !
703   END SUBROUTINE ts_rst
704
705#else
706   !!----------------------------------------------------------------------
707   !!   Default case :   Empty module   No standart free surface cst volume
708   !!----------------------------------------------------------------------
709CONTAINS
710   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function
711      dyn_spg_ts_alloc = 0
712   END FUNCTION dyn_spg_ts_alloc
713   SUBROUTINE dyn_spg_ts( kt )            ! Empty routine
714      INTEGER, INTENT(in) :: kt
715      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
716   END SUBROUTINE dyn_spg_ts
717   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine
718      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
719      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
720      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
721   END SUBROUTINE ts_rst   
722#endif
723   
724   !!======================================================================
725END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.