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

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

source: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 363

Last change on this file since 363 was 363, checked in by opalod, 18 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 24.3 KB
Line 
1MODULE dynspg_ts
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_ts  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if ( defined key_dynspg_ts && ! defined key_autotasking ) ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_ts'     free surface cst volume with time splitting
9   !!   NOT 'key_autotasking'                      k-j-i loop (vector opt.)
10   !!----------------------------------------------------------------------
11   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
12   !!                 splitting scheme and add to the general trend
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! physical constants
18   USE ocesbc          ! ocean surface boundary condition
19   USE dynvor          ! vorticity term
20   USE obc_oce         ! Lateral open boundary condition
21   USE lib_mpp         ! distributed memory computing library
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl          ! Print control
24   USE in_out_manager  ! I/O manager
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Accessibility
30   PUBLIC dyn_spg_ts  ! routine called by step.F90
31
32   !! * Module variables
33      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop
34         sshn_b, sshb_b,               &  ! sea surface heigth (now, before)
35         un_b  , vn_b                     ! vertically integrated horizontal velocities (now)
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !!  OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Header$
43   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE dyn_spg_ts( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  routine dyn_spg_ts  ***
51      !!
52      !! ** Purpose :   Compute the now trend due to the surface pressure
53      !!      gradient in case of free surface formulation with time-splitting.
54      !!      Add it to the general trend of momentum equation.
55      !!      Compute the free surface.
56      !!
57      !! ** Method  :   Free surface formulation with time-splitting
58      !!      -1- Save the vertically integrated trend. This general trend is
59      !!          held constant over the barotropic integration.
60      !!          The Coriolis force is removed from the general trend as the
61      !!          surface gradient and the Coriolis force are updated within
62      !!          the barotropic integration.
63      !!      -2- Barotropic loop : updates of sea surface height (zssha_e) and
64      !!          barotropic transports (zua_e and zva_e) through barotropic
65      !!          momentum and continuity integration. Barotropic former
66      !!          variables are time averaging over the full barotropic cycle
67      !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b
68      !!          and zvX_b (X specifying after, now or before).
69      !!      -3- Update of sea surface height from time averaged barotropic
70      !!          variables.
71      !!        - apply lateral boundary conditions on sshn.
72      !!      -4- The new general trend becomes :
73      !!          ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H
74      !!
75      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
76      !!
77      !! References :
78      !!   Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
79      !!
80      !! History :
81      !!   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code
82      !!        !  05-11  (V. Garnier, G. Madec)  optimization
83      !!---------------------------------------------------------------------
84      !! * Arguments
85      INTEGER, INTENT( in )  ::   kt           ! ocean time-step index
86
87      !! * Local declarations
88      INTEGER  ::  ji, jj, jk, jit             ! dummy loop indices
89      INTEGER  ::  icycle                      ! temporary scalar
90      REAL(wp) ::                           &
91         zraur, zcoef, z2dt_e, z2dt_b, zfac25,   &  ! temporary scalars
92         zfact1, zspgu, zcubt, zx1, zy1,    &  !     "        "
93         zfact2, zspgv, zcvbt, zx2, zy2        !     "        "
94      REAL(wp), DIMENSION(jpi,jpj) ::       &
95         zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays
96         zua, zva, zub, zvb,                &  !     "        "
97         zssha_b, zua_b, zva_b,             &  !     "        "
98         zsshb_e, zub_e, zvb_e,             &  !     "        "
99         zsshn_e, zun_e, zvn_e,             &  !     "        "
100         zssha_e, zua_e, zva_e                 !     "        "
101      REAL(wp), DIMENSION(jpi,jpj),SAVE ::  &
102         ztnw, ztne, ztsw, ztse
103      !!----------------------------------------------------------------------
104
105      ! Arrays initialization
106      ! ---------------------
107      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0   ;   zua_e(:,:) = 0.e0
108      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0   ;   zva_e(:,:) = 0.e0
109      zhdiv(:,:) = 0.e0
110
111
112      IF( kt == nit000 ) THEN
113
114         IF(lwp) WRITE(numout,*)
115         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
116         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
117         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2*rdt/rdtbt )
118
119         IF( ln_rstart ) THEN
120            ! initialize barotropic specific arrays
121            sshb_b(:,:) = sshb(:,:)
122            sshn_b(:,:) = sshn(:,:)
123            un_b(:,:)   = 0.e0
124            vn_b(:,:)   = 0.e0
125            ! vertical sum
126            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
127               DO jk = 1, jpkm1
128                  DO ji = 1, jpij
129                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
130                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
131                  END DO
132               END DO
133            ELSE                             ! No  vector opt.
134               DO jk = 1, jpkm1
135                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
136                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
137               END DO
138            ENDIF
139         ENDIF
140         zssha_e(:,:) = sshn(:,:)
141         zua_e  (:,:) = un_b(:,:)
142         zva_e  (:,:) = vn_b(:,:)
143
144         IF( ln_dynvor_een ) THEN
145            ztne(1,:) = 0.e0   ;   ztnw(1,:) = 0.e0   ;   ztse(1,:) = 0.e0   ;   ztsw(1,:) = 0.e0
146            DO jj = 2, jpj
147               DO ji = fs_2, jpi   ! vector opt.
148                  ztne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
149                  ztnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
150                  ztse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
151                  ztsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
152               END DO
153            END DO
154         ENDIF
155
156      ENDIF
157   
158      ! Local constant initialization
159      ! --------------------------------
160      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
161      IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
162      zfact1 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
163      zfact2 = 0.5 * 0.5
164      zraur  = 1. / rauw                                   ! 1 / volumic mass of pure water
165     
166      ! -----------------------------------------------------------------------------
167      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
168      ! -----------------------------------------------------------------------------
169
170      ! Vertically integrated quantities
171      ! --------------------------------
172      zua(:,:) = 0.e0
173      zva(:,:) = 0.e0
174      zub(:,:) = 0.e0
175      zvb(:,:) = 0.e0
176      zwx(:,:) = 0.e0
177      zwy(:,:) = 0.e0
178
179      ! vertical sum
180      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
181         DO jk = 1, jpkm1
182            DO ji = 1, jpij
183               !                                                           ! Vertically integrated momentum trends
184               zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk)
185               zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
186               !                                                           ! Vertically integrated transports (before)
187               zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk)
188               zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk)
189               !                                                           ! Planetary vorticity transport fluxes (now)
190               zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
191               zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk)
192            END DO
193         END DO
194      ELSE                             ! No  vector opt.
195         DO jk = 1, jpkm1
196            !                                                           ! Vertically integrated momentum trends
197            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
198            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
199            !                                                           ! Vertically integrated transports (before)
200            zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk)
201            zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk)
202            !                                                           ! Planetary vorticity (now)
203            zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
204            zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
205         END DO
206      ENDIF
207
208      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
209         DO jj = 2, jpjm1
210            DO ji = fs_2, fs_jpim1   ! vector opt.
211               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
212               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
213               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
214               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
215               ! energy conserving formulation for planetary vorticity term
216               zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
217               zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
218            END DO
219         END DO
220
221      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
222         DO jj = 2, jpjm1
223            DO ji = fs_2, fs_jpim1   ! vector opt.
224               zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
225                              + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
226               zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
227                              + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
228               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
229               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
230            END DO
231         END DO
232
233      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
234         zfac25 = 0.25
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1   ! vector opt.
237               zcu(ji,jj) = + zfac25 / e1u(ji,jj)   &
238                  &       * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
239                  &          + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
240               zcv(ji,jj) = - zfac25 / e2v(ji,jj)   &
241                  &       * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
242                  &          + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
243            END DO
244         END DO
245
246      ENDIF
247
248
249      ! Remove barotropic trend from general momentum trend
250      ! ---------------------------------------------------
251      DO jk = 1 , jpkm1
252         DO jj = 2, jpjm1
253            DO ji = fs_2, fs_jpim1   ! vector opt.
254               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
255               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
256            END DO
257         END DO
258      END DO
259
260      ! Remove coriolis term from barotropic trend
261      ! ------------------------------------------
262      DO jj = 2, jpjm1
263         DO ji = fs_2, fs_jpim1
264            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
265            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
266         END DO
267      END DO
268
269      ! -----------------------------------------------------------------------
270      !  Phase 2 : Integration of the barotropic equations with time splitting
271      ! -----------------------------------------------------------------------
272
273      ! Initialisations
274      !----------------
275      ! Number of iteration of the barotropic loop
276      icycle = FLOOR( z2dt_b / rdtbt )
277
278      ! variables for the barotropic equations
279      zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
280      zsshn_e(:,:) = sshn_b(:,:)
281      zub_e(:,:) = un_b(:,:)       ! barotropic transports issued from the barotropic equations (before and now)
282      zvb_e(:,:) = vn_b(:,:)
283      zun_e(:,:) = un_b(:,:)
284      zvn_e(:,:) = vn_b(:,:)
285      zssha_b(:,:) = sshn(:,:)        ! time averaged variables over all sub-timesteps
286      zua_b(:,:) = un_b(:,:)   
287      zva_b(:,:) = vn_b(:,:)
288
289      ! Barotropic integration over 2 baroclinic time steps
290      ! ---------------------------------------------------
291
292      !                                                    ! ==================== !
293      DO jit = 1, icycle                                   !  sub-time-step loop  !
294         !                                                 ! ==================== !
295
296         z2dt_e = 2. * rdtbt
297         IF ( jit == 1 )   z2dt_e = rdtbt
298
299         ! Horizontal divergence of barotropic transports
300         !--------------------------------------------------
301         DO jj = 2, jpjm1
302            DO ji = fs_2, fs_jpim1   ! vector opt.
303               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              &
304                  &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              &
305                  &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              &
306                  &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          &
307                  &           / (e1t(ji,jj)*e2t(ji,jj))
308            END DO
309         END DO
310
311#if defined key_obc
312         ! open boundaries (div must be zero behind the open boundary)
313         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
314         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
315         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
316         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
317         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
318#endif
319
320         ! Sea surface height from the barotropic system
321         !----------------------------------------------
322         DO jj = 2, jpjm1
323            DO ji = fs_2, fs_jpim1   ! vector opt.
324               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
325            &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
326            END DO
327         END DO
328
329         ! evolution of the barotropic transport ( following the vorticity scheme used)
330         ! ----------------------------------------------------------------------------
331         zwx(:,:) = e2u(:,:) * zun_e(:,:)
332         zwy(:,:) = e1v(:,:) * zvn_e(:,:)
333
334         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
335            DO jj = 2, jpjm1
336               DO ji = fs_2, fs_jpim1   ! vector opt.
337                  ! surface pressure gradient
338                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
339                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
340                  ! energy conserving formulation for planetary vorticity term
341                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
342                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
343                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
344                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
345                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
346                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
347                  ! after transports
348                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
349                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
350               END DO
351            END DO
352
353         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
354            DO jj = 2, jpjm1
355               DO ji = fs_2, fs_jpim1   ! vector opt.
356                  ! surface pressure gradient
357                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
358                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
359                  ! enstrophy conserving formulation for planetary vorticity term
360                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
361                                 + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
362                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
363                                 + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
364                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
365                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
366                  ! after transports
367                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
368                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
369               END DO
370            END DO
371
372         ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
373            zfac25 = 0.25
374            DO jj = 2, jpjm1
375               DO ji = fs_2, fs_jpim1   ! vector opt.
376                  ! surface pressure gradient
377                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
378                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
379                  ! energy/enstrophy conserving formulation for planetary vorticity term
380                  zcubt = + zfac25 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
381                     &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
382                  zcvbt = - zfac25 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
383                     &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
384                  ! after transports
385                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
386                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
387               END DO
388            END DO
389         ENDIF
390
391         ! ... Boundary conditions on zua_e, zva_e, zssha_e
392         CALL lbc_lnk( zua_e, 'U', -1. )
393         CALL lbc_lnk( zva_e, 'V', -1. )
394         CALL lbc_lnk( zssha_e, 'T', 1. )
395
396         ! temporal sum
397         !-------------
398         zssha_b(:,:) = zssha_b(:,:) + zssha_e(:,:)
399         zua_b  (:,:) = zua_b  (:,:) + zua_e  (:,:)
400         zva_b  (:,:) = zva_b  (:,:) + zva_e  (:,:) 
401
402         ! Time filter and swap of dynamics arrays
403         ! ---------------------------------------
404         IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
405            zsshb_e(:,:) = zsshn_e(:,:)
406            zub_e  (:,:) = zun_e  (:,:)
407            zvb_e  (:,:) = zvn_e  (:,:)
408            zsshn_e(:,:) = zssha_e(:,:)
409            zun_e  (:,:) = zua_e  (:,:)
410            zvn_e  (:,:) = zva_e  (:,:)
411         ELSE                                        ! Asselin filtering
412            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * zsshn_e(:,:)
413            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * zun_e  (:,:)
414            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * zvn_e  (:,:)
415            zsshn_e(:,:) = zssha_e(:,:)
416            zun_e  (:,:) = zua_e  (:,:)
417            zvn_e  (:,:) = zva_e  (:,:)
418         ENDIF
419
420         !                                                 ! ==================== !
421      END DO                                               !        end loop      !
422      !                                                    ! ==================== !
423
424
425      ! Time average of after barotropic variables
426      zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
427      zssha_b(:,:) = zcoef * zssha_b(:,:) 
428      zua_b  (:,:) = zcoef * zua_b  (:,:) 
429      zva_b  (:,:) = zcoef * zva_b  (:,:) 
430     
431
432      ! ---------------------------------------------------------------------------
433      ! Phase 3 : Update sea surface height from time averaged barotropic variables
434      ! ---------------------------------------------------------------------------
435
436 
437      ! Horizontal divergence of time averaged barotropic transports
438      !-------------------------------------------------------------
439      DO jj = 2, jpjm1
440         DO ji = fs_2, fs_jpim1   ! vector opt.
441            zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
442           &                +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
443           &             / ( e1t(ji,jj) * e2t(ji,jj) )
444         END DO
445      END DO
446
447#if defined key_obc
448      ! open boundaries (div must be zero behind the open boundary)
449      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
450      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
451      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
452      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
453      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
454#endif
455
456      ! sea surface height
457      !-------------------
458      sshb(:,:) = sshn(:,:)
459      sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
460
461      ! ... Boundary conditions on sshn
462      CALL lbc_lnk( sshn, 'T', 1. )
463
464
465      ! -----------------------------------------------------------------------------
466      ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
467      ! -----------------------------------------------------------------------------
468
469      ! Swap on time averaged barotropic variables
470      ! ------------------------------------------
471      sshb_b(:,:) = sshn_b (:,:)
472      sshn_b(:,:) = zssha_b(:,:)
473      un_b  (:,:) = zua_b  (:,:) 
474      vn_b  (:,:) = zva_b  (:,:) 
475   
476      ! add time averaged barotropic coriolis and surface pressure gradient
477      ! terms to the general momentum trend
478      ! --------------------------------------------------------------------
479      DO jk=1,jpkm1
480         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b
481         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b
482      END DO
483
484      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
485         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
486      ENDIF
487     
488   END SUBROUTINE dyn_spg_ts
489#else
490   !!----------------------------------------------------------------------
491   !!   Default case :   Empty module   No standart free surface cst volume
492   !!----------------------------------------------------------------------
493CONTAINS
494   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
495      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
496   END SUBROUTINE dyn_spg_ts
497#endif
498   
499   !!======================================================================
500END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.