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

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

nemo_v1_update_035 : CT : OBCs adapted to the 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: 25.1 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 obcdta          ! open boundary condition data     
20   USE obcfla          ! Flather open boundary condition 
21   USE dynvor          ! vorticity term
22   USE obc_oce         ! Lateral open boundary condition
23   USE lib_mpp         ! distributed memory computing library
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE prtctl          ! Print control
26   USE dynspg_oce      ! surface pressure gradient variables
27   USE in_out_manager  ! I/O manager
28
29   IMPLICIT NONE
30   PRIVATE
31
32   !! * Accessibility
33   PUBLIC dyn_spg_ts  ! routine called by step.F90
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !!  OPA 9.0 , LOCEAN-IPSL (2005)
40   !! $Header$
41   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
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 (ssha_e) and
62      !!          barotropic transports (ua_e and va_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         zun_e, zvn_e                          !     "        "
98      REAL(wp), DIMENSION(jpi,jpj),SAVE ::  &
99         ztnw, ztne, ztsw, ztse
100      !!----------------------------------------------------------------------
101
102      ! Arrays initialization
103      ! ---------------------
104      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0
105      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0
106      zhdiv(:,:) = 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         ssha_e(:,:) = sshn(:,:)
138         ua_e(:,:)   = un_b(:,:)
139         va_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      sshn_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      ! set ssh corrections to 0
287      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
288#if defined key_obc
289      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
290      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
291      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
292      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
293#endif
294
295      ! Barotropic integration over 2 baroclinic time steps
296      ! ---------------------------------------------------
297
298      !                                                    ! ==================== !
299      DO jit = 1, icycle                                   !  sub-time-step loop  !
300         !                                                 ! ==================== !
301
302         z2dt_e = 2. * rdtbt
303         IF ( jit == 1 )   z2dt_e = rdtbt
304
305         ! Time interpolation of open boundary condition data
306         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
307
308         ! Horizontal divergence of barotropic transports
309         !--------------------------------------------------
310         DO jj = 2, jpjm1
311            DO ji = fs_2, fs_jpim1   ! vector opt.
312               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              &
313                  &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              &
314                  &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              &
315                  &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          &
316                  &           / (e1t(ji,jj)*e2t(ji,jj))
317            END DO
318         END DO
319
320#if defined key_obc
321         ! open boundaries (div must be zero behind the open boundary)
322         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
323         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0      ! east
324         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0      ! west
325         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
326         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0      ! south
327#endif
328
329         ! Sea surface height from the barotropic system
330         !----------------------------------------------
331         DO jj = 2, jpjm1
332            DO ji = fs_2, fs_jpim1   ! vector opt.
333               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
334            &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
335            END DO
336         END DO
337
338         ! evolution of the barotropic transport ( following the vorticity scheme used)
339         ! ----------------------------------------------------------------------------
340         zwx(:,:) = e2u(:,:) * zun_e(:,:)
341         zwy(:,:) = e1v(:,:) * zvn_e(:,:)
342
343         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
344            DO jj = 2, jpjm1
345               DO ji = fs_2, fs_jpim1   ! vector opt.
346                  ! surface pressure gradient
347                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
348                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
349                  ! energy conserving formulation for planetary vorticity term
350                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
351                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
352                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
353                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
354                  zcubt = zfact2 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 )
355                  zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
356                  ! after transports
357                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
358                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
359               END DO
360            END DO
361
362         ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
363            DO jj = 2, jpjm1
364               DO ji = fs_2, fs_jpim1   ! vector opt.
365                  ! surface pressure gradient
366                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
367                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
368                  ! enstrophy conserving formulation for planetary vorticity term
369                  zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
370                                 + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
371                  zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
372                                 + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
373                  zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
374                  zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
375                  ! after transports
376                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
377                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
378               END DO
379            END DO
380
381         ELSEIF ( ln_dynvor_een ) THEN                    ! energy and enstrophy conserving scheme
382            zfac25 = 0.25
383            DO jj = 2, jpjm1
384               DO ji = fs_2, fs_jpim1   ! vector opt.
385                  ! surface pressure gradient
386                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
387                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
388                  ! energy/enstrophy conserving formulation for planetary vorticity term
389                  zcubt = + zfac25 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
390                     &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
391                  zcvbt = - zfac25 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
392                     &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
393                  ! after transports
394                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
395                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
396               END DO
397            END DO
398         ENDIF
399
400         ! Flather's boundary condition for the barotropic loop :
401         !         - Update sea surface height on each open boundary
402         !         - Correct the barotropic transports
403         IF( lk_obc )   CALL obc_fla_ts( kt )
404
405
406         ! ... Boundary conditions on ua_e, va_e, ssha_e
407         CALL lbc_lnk( ua_e  , 'U', -1. )
408         CALL lbc_lnk( va_e  , 'V', -1. )
409         CALL lbc_lnk( ssha_e, 'T',  1. )
410
411         ! temporal sum
412         !-------------
413         zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:)
414         zua_b  (:,:) = zua_b  (:,:) + ua_e  (:,:)
415         zva_b  (:,:) = zva_b  (:,:) + va_e  (:,:) 
416
417         ! Time filter and swap of dynamics arrays
418         ! ---------------------------------------
419         IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
420            zsshb_e(:,:) = sshn_e(:,:)
421            zub_e  (:,:) = zun_e (:,:)
422            zvb_e  (:,:) = zvn_e (:,:)
423            sshn_e (:,:) = ssha_e(:,:)
424            zun_e  (:,:) = ua_e  (:,:)
425            zvn_e  (:,:) = va_e  (:,:)
426         ELSE                                        ! Asselin filtering
427            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
428            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:)
429            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:)
430            sshn_e (:,:) = ssha_e(:,:)
431            zun_e  (:,:) = ua_e  (:,:)
432            zvn_e  (:,:) = va_e  (:,:)
433         ENDIF
434
435         !                                                 ! ==================== !
436      END DO                                               !        end loop      !
437      !                                                    ! ==================== !
438
439
440      ! Time average of after barotropic variables
441      zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
442      zssha_b(:,:) = zcoef * zssha_b(:,:) 
443      zua_b  (:,:) = zcoef *  zua_b (:,:) 
444      zva_b  (:,:) = zcoef *  zva_b (:,:) 
445#if defined key_obc
446         IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
447         IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
448         IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
449         IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
450#endif
451     
452
453      ! ---------------------------------------------------------------------------
454      ! Phase 3 : Update sea surface height from time averaged barotropic variables
455      ! ---------------------------------------------------------------------------
456
457 
458      ! Horizontal divergence of time averaged barotropic transports
459      !-------------------------------------------------------------
460      DO jj = 2, jpjm1
461         DO ji = fs_2, fs_jpim1   ! vector opt.
462            zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
463           &                +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
464           &             / ( e1t(ji,jj) * e2t(ji,jj) )
465         END DO
466      END DO
467
468#if defined key_obc
469      ! open boundaries (div must be zero behind the open boundary)
470      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
471      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
472      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
473      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
474      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
475#endif
476
477      ! sea surface height
478      !-------------------
479      sshb(:,:) = sshn(:,:)
480      sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
481
482      ! ... Boundary conditions on sshn
483      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
484
485
486      ! -----------------------------------------------------------------------------
487      ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
488      ! -----------------------------------------------------------------------------
489
490      ! Swap on time averaged barotropic variables
491      ! ------------------------------------------
492      sshb_b(:,:) = sshn_b (:,:)
493      sshn_b(:,:) = zssha_b(:,:)
494      un_b  (:,:) = zua_b  (:,:) 
495      vn_b  (:,:) = zva_b  (:,:) 
496   
497      ! add time averaged barotropic coriolis and surface pressure gradient
498      ! terms to the general momentum trend
499      ! --------------------------------------------------------------------
500      DO jk=1,jpkm1
501         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b
502         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b
503      END DO
504
505      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
506         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
507      ENDIF
508     
509   END SUBROUTINE dyn_spg_ts
510#else
511   !!----------------------------------------------------------------------
512   !!   Default case :   Empty module   No standart free surface cst volume
513   !!----------------------------------------------------------------------
514CONTAINS
515   SUBROUTINE dyn_spg_ts( kt )       ! Empty routine
516      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt
517   END SUBROUTINE dyn_spg_ts
518#endif
519   
520   !!======================================================================
521END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.