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

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