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

source: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4376

Last change on this file since 4376 was 4376, checked in by hliu, 10 years ago

updating wetting/drying filter for explicit time_splitting

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