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

Last change on this file since 455 was 455, checked in by opalod, 15 years ago

nemo_v1_update_048:RB: reorganization of dynamics part, in addition change atsk to jki, suppress dynhpg_atsk.F90 dynzdf_imp_atsk.F90 dynzdf_iso.F90

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