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

Last change on this file since 359 was 359, 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
110
111      IF( kt == nit000 ) THEN
112
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
115         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
116         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2*rdt/rdtbt )
117
118         IF( ln_rstart ) THEN
119            ! initialize barotropic specific arrays
120            sshb_b(:,:) = sshb(:,:)
121            sshn_b(:,:) = sshn(:,:)
122            un_b(:,:)   = 0.e0
123            vn_b(:,:)   = 0.e0
124            ! vertical sum
125            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
126               DO jk = 1, jpkm1
127                  DO ji = 1, jpij
128                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
129                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
130                  END DO
131               END DO
132            ELSE                             ! No  vector opt.
133               DO jk = 1, jpkm1
134                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
135                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
136               END DO
137            ENDIF
138         ENDIF
139         zssha_e(:,:) = sshn(:,:)
140         zua_e  (:,:) = un_b(:,:)
141         zva_e  (:,:) = vn_b(:,:)
142
143         IF( ln_dynvor_een ) THEN
144            ztne(1,:) = 0.e0   ;   ztnw(1,:) = 0.e0   ;   ztse(1,:) = 0.e0   ;   ztsw(1,:) = 0.e0
145            DO jj = 2, jpj
146               DO ji = fs_2, jpi   ! vector opt.
147                  ztne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
148                  ztnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
149                  ztse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3.
150                  ztsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3.
151               END DO
152            END DO
153         ENDIF
154
155      ENDIF
156   
157      ! Local constant initialization
158      ! --------------------------------
159      z2dt_b = 2.0 * rdt                                    ! baroclinic time step
160      IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
161      zfact1 = 0.5 * 0.25                                     ! coefficient for vorticity estimates
162      zfact2 = 0.5 * 0.5
163      zraur  = 1. / rauw                                   ! 1 / volumic mass of pure water
164     
165      ! -----------------------------------------------------------------------------
166      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
167      ! -----------------------------------------------------------------------------
168
169      ! Vertically integrated quantities
170      ! --------------------------------
171      zua(:,:) = 0.e0
172      zva(:,:) = 0.e0
173      zub(:,:) = 0.e0
174      zvb(:,:) = 0.e0
175      zwx(:,:) = 0.e0
176      zwy(:,:) = 0.e0
177
178      ! vertical sum
179      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
180         DO jk = 1, jpkm1
181            DO ji = 1, jpij
182               !                                                           ! Vertically integrated momentum trends
183               zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk)
184               zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
185               !                                                           ! Vertically integrated transports (before)
186               zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk)
187               zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk)
188               !                                                           ! Planetary vorticity transport fluxes (now)
189               zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
190               zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk)
191            END DO
192         END DO
193      ELSE                             ! No  vector opt.
194         DO jk = 1, jpkm1
195            !                                                           ! Vertically integrated momentum trends
196            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
197            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
198            !                                                           ! Vertically integrated transports (before)
199            zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk)
200            zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk)
201            !                                                           ! Planetary vorticity (now)
202            zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
203            zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
204         END DO
205      ENDIF
206
207      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
208         DO jj = 2, jpjm1
209            DO ji = fs_2, fs_jpim1   ! vector opt.
210               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
211               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
212               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
213               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
214               ! energy conserving formulation for planetary vorticity term
215               zcu(ji,jj) = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
216               zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
217            END DO
218         END DO
219
220      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
221         DO jj = 2, jpjm1
222            DO ji = fs_2, fs_jpim1   ! vector opt.
223               zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
224                              + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
225               zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
226                              + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
227               zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
228               zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
229            END DO
230         END DO
231
232      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
233         zfac25 = 0.25
234         DO jj = 2, jpjm1
235            DO ji = fs_2, fs_jpim1   ! vector opt.
236               zcu(ji,jj) = + zfac25 / e1u(ji,jj)   &
237                  &       * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
238                  &          + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
239               zcv(ji,jj) = - zfac25 / e2v(ji,jj)   &
240                  &       * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
241                  &          + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
242            END DO
243         END DO
244
245      ENDIF
246
247
248      ! Remove barotropic trend from general momentum trend
249      ! ---------------------------------------------------
250      DO jk = 1 , jpkm1
251         DO jj = 2, jpjm1
252            DO ji = fs_2, fs_jpim1   ! vector opt.
253               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
254               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
255            END DO
256         END DO
257      END DO
258
259      ! Remove coriolis term from barotropic trend
260      ! ------------------------------------------
261      DO jj = 2, jpjm1
262         DO ji = fs_2, fs_jpim1
263            zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
264            zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
265         END DO
266      END DO
267
268      ! -----------------------------------------------------------------------
269      !  Phase 2 : Integration of the barotropic equations with time splitting
270      ! -----------------------------------------------------------------------
271
272      ! Initialisations
273      !----------------
274      ! Number of iteration of the barotropic loop
275      icycle = FLOOR( z2dt_b / rdtbt )
276
277      ! variables for the barotropic equations
278      zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
279      zsshn_e(:,:) = sshn_b(:,:)
280      zub_e(:,:) = un_b(:,:)       ! barotropic transports issued from the barotropic equations (before and now)
281      zvb_e(:,:) = vn_b(:,:)
282      zun_e(:,:) = un_b(:,:)
283      zvn_e(:,:) = vn_b(:,:)
284      zssha_b(:,:) = sshn(:,:)        ! time averaged variables over all sub-timesteps
285      zua_b(:,:) = un_b(:,:)   
286      zva_b(:,:) = vn_b(:,:)
287
288      ! Barotropic integration over 2 baroclinic time steps
289      ! ---------------------------------------------------
290
291      !                                                    ! ==================== !
292      DO jit = 1, icycle                                   !  sub-time-step loop  !
293         !                                                 ! ==================== !
294
295         z2dt_e = 2. * rdtbt
296         IF ( jit == 1 )   z2dt_e = rdtbt
297
298         ! Horizontal divergence of barotropic transports
299         !--------------------------------------------------
300         DO jj = 2, jpjm1
301            DO ji = fs_2, fs_jpim1   ! vector opt.
302               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              &
303                  &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              &
304                  &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              &
305                  &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          &
306                  &           / (e1t(ji,jj)*e2t(ji,jj))
307            END DO
308         END DO
309
310#if defined key_obc
311         ! open boundaries (div must be zero behind the open boundary)
312         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
313         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
314         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
315         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
316         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
317#endif
318
319         ! Sea surface height from the barotropic system
320         !----------------------------------------------
321         DO jj = 2, jpjm1
322            DO ji = fs_2, fs_jpim1   ! vector opt.
323               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
324            &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
325            END DO
326         END DO
327
328         ! evolution of the barotropic transport ( following the vorticity scheme used)
329         ! ----------------------------------------------------------------------------
330         zwx(:,:) = e2u(:,:) * zun_e(:,:)
331         zwy(:,:) = e1v(:,:) * zvn_e(:,:)
332
333         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
334            DO jj = 2, jpjm1
335               DO ji = fs_2, fs_jpim1   ! vector opt.
336                  ! surface pressure gradient
337                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
338                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
339                  ! energy conserving formulation for planetary vorticity term
340                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
341                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
342                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
343                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
344                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
345                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
346                  ! after transports
347                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
348                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
349               END DO
350            END DO
351
352         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
353            DO jj = 2, jpjm1
354               DO ji = fs_2, fs_jpim1   ! vector opt.
355                  ! surface pressure gradient
356                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
357                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
358                  ! enstrophy conserving formulation for planetary vorticity term
359                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
360                                 + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
361                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
362                                 + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
363                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
364                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
365                  ! after transports
366                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
367                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
368               END DO
369            END DO
370
371         ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
372            zfac25 = 0.25
373            DO jj = 2, jpjm1
374               DO ji = fs_2, fs_jpim1   ! vector opt.
375                  ! surface pressure gradient
376                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
377                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
378                  ! energy/enstrophy conserving formulation for planetary vorticity term
379                  zcubt = + zfac25 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
380                     &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
381                  zcvbt = - zfac25 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
382                     &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
383                  ! after transports
384                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
385                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
386               END DO
387            END DO
388         ENDIF
389
390         ! ... Boundary conditions on zua_e, zva_e, zssha_e
391         CALL lbc_lnk( zua_e, 'U', -1. )
392         CALL lbc_lnk( zva_e, 'V', -1. )
393         CALL lbc_lnk( zssha_e, 'T', 1. )
394
395         ! temporal sum
396         !-------------
397         zssha_b(:,:) = zssha_b(:,:) + zssha_e(:,:)
398         zua_b  (:,:) = zua_b  (:,:) + zua_e  (:,:)
399         zva_b  (:,:) = zva_b  (:,:) + zva_e  (:,:) 
400
401         ! Time filter and swap of dynamics arrays
402         ! ---------------------------------------
403         IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
404            zsshb_e(:,:) = zsshn_e(:,:)
405            zub_e  (:,:) = zun_e  (:,:)
406            zvb_e  (:,:) = zvn_e  (:,:)
407            zsshn_e(:,:) = zssha_e(:,:)
408            zun_e  (:,:) = zua_e  (:,:)
409            zvn_e  (:,:) = zva_e  (:,:)
410         ELSE                                        ! Asselin filtering
411            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * zsshn_e(:,:)
412            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * zun_e  (:,:)
413            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * zvn_e  (:,:)
414            zsshn_e(:,:) = zssha_e(:,:)
415            zun_e  (:,:) = zua_e  (:,:)
416            zvn_e  (:,:) = zva_e  (:,:)
417         ENDIF
418
419         !                                                 ! ==================== !
420      END DO                                               !        end loop      !
421      !                                                    ! ==================== !
422
423
424      ! Time average of after barotropic variables
425      zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
426      zssha_b(:,:) = zcoef * zssha_b(:,:) 
427      zua_b  (:,:) = zcoef * zua_b  (:,:) 
428      zva_b  (:,:) = zcoef * zva_b  (:,:) 
429     
430
431      ! ---------------------------------------------------------------------------
432      ! Phase 3 : Update sea surface height from time averaged barotropic variables
433      ! ---------------------------------------------------------------------------
434
435 
436      ! Horizontal divergence of time averaged barotropic transports
437      !-------------------------------------------------------------
438      DO jj = 2, jpjm1
439         DO ji = fs_2, fs_jpim1   ! vector opt.
440            zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
441           &                +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
442           &             / ( e1t(ji,jj) * e2t(ji,jj) )
443         END DO
444      END DO
445
446#if defined key_obc
447      ! open boundaries (div must be zero behind the open boundary)
448      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
449      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
450      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
451      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
452      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
453#endif
454
455      ! sea surface height
456      !-------------------
457      sshb(:,:) = sshn(:,:)
458      sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
459
460      ! ... Boundary conditions on sshn
461      CALL lbc_lnk( sshn, 'T', 1. )
462
463
464      ! -----------------------------------------------------------------------------
465      ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
466      ! -----------------------------------------------------------------------------
467
468      ! Swap on time averaged barotropic variables
469      ! ------------------------------------------
470      sshb_b(:,:) = sshn_b (:,:)
471      sshn_b(:,:) = zssha_b(:,:)
472      un_b  (:,:) = zua_b  (:,:) 
473      vn_b  (:,:) = zva_b  (:,:) 
474   
475      ! add time averaged barotropic coriolis and surface pressure gradient
476      ! terms to the general momentum trend
477      ! --------------------------------------------------------------------
478      DO jk=1,jpkm1
479         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b
480         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b
481      END DO
482
483      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
484         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
485      ENDIF
486     
487   END SUBROUTINE dyn_spg_ts
488#else
489   !!----------------------------------------------------------------------
490   !!   Default case :   Empty module   No standart free surface cst volume
491   !!----------------------------------------------------------------------
492CONTAINS
493   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
494      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
495   END SUBROUTINE dyn_spg_ts
496#endif
497   
498   !!======================================================================
499END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.